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ABSTRACT 


Two  computer  programs  for  simulating  countercurrent 
distribution  are  presented.     The  programs  can  simulate  frontal 
and  elution  modes  of  operation.     Phase  volumes  and  composi- 
tions in  all  tubes  at  any  point  in  the  distribution  may  be  printed 
out.     Effluent  information  is  printed  at  the  end  of  a  run.     Vari- 
able partition  coefficients,   two-phase  flow,   and  many  other  non- 
ideal  phenomena  can  be  simulated.     Program  I  can  be  used  for 
a  quick  approximation  and  is  accurate  in  dilute  and  near  ideal 
situations.     Program  II,    incorporating  phase  diagrams,   gives 
good  agreement  with  experimental  evidence  but  is  more  diffi- 
cult to  set  up. 
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COMPUTER    PROGRAMS   FOR   COUNTERCURRENT   DISTRIBUTION 


V.   G.    Martin,    R.    A.    Barford,    C.    R.    Eddy,    and  H.    L.    Rothbart 


A  recent  publication*  considered  nonideal  phenomena  encountered  in  counter- 
current  distribution.     We  found  it  necessary  to  simulate  these  processes  on  the  com- 
puter.    Two  programs  (I  and  II)  were  written  in  FORTRAN  IV  for   an  IBM    1130** 
equipped  with  card  reader-punch,    disk,  8  K  core. 


PROGRAM   I 

Program  I  is  less  fundamental  than  Program  II,    but  it  is  a  good  approxima- 
tion under  the  following  conditions:    The  solutes  effectively  have  no  volume;  there  is 
no  solute-solute  interaction;  the  partition  coefficients  are  constant  or  vary  as  a  func- 
tion of  total  concentration;  solvent  volumes  are  constant.     These  assumptions  are 
reasonable  in  dilute  solutions. 

Two  simulations  may  be  done  simultaneously.     The  partition  coefficient  and 
concentrations  for  the  first  solute  are  designated  by  variables  in  the. program  begin- 
ning with  the  letter  C.     The  variables  associated  with  the  second  solute  begin  with  the 
letter  X.     An  apparatus  of  1  to  200  tubes  may  be  simulated.     The  number  of  tubes  in 
the  apparatus  and  the  volume  information  for  both  systems  are  the  same. 

For  each  system,    one  or  more  inputs  of  solute  may  be  desired.     An  elution 
profile  results  from  a  single  input  of  solute.     If  enough  increments  of  solute  are  add- 
ed to  tube  zero,   the  maximum  equilibrium  concentration  will  be  reached  in  the  last 
tube  and  a  frontal  output  profile  will  result.    A  number  of  inputs  between  two  and  the 
number  required  to  just  reach  the  maximum  equilibrium  concentration  will  give  an 
output  profile  somewhere  between  an  elution  and  a  frontal.     This  region  may  be  inter- 
esting in  certain  systems,    and  it  is  particularly  suited  to  computer  simulation  since 
this  intermediate  region  is  difficult  to  describe  mathematically. 

Solute  may  be  introduced  to  the  instrument  by  two  different  procedures.     If 
data  switch  '14'  is  on,    the  solute  is  placed  in  one  or  more  tubes  before  the  simulation 
is  started.     This  batch  loading  reduces  the  number  of  effective  tubes  in  the  instru- 


*Rothbart,    H.    L.  ,    Barford,    R.    A.,    Martin,    V.    G.  ,    Bertsch,    R.    J.,    and 
Eddy,   C.   R.     Nonideal  Phenomena  in  Countercurrent  Distribution.     Separation  Sci- 
ence  [In  press.  ] 

**Trade  names  are  used  in  this  publication  solely  for  the  purpose  of  providing 
specific  information.     Mention  of  a  trade  name  does  not  constitute  a  guarantee  or 
warranty  of  the  product  by  the  U.   S.    Department  of  Agriculture  or  an  endorsement 
by  the  Department  over  other  products  not  mentioned. 


ment.     If  data  switch  '14'  is  off,   the  solute  is  added  to  the  zeroth  tube  in  one  or  more 
increments  with  a  transfer  occurring  between  each  addition.     This  sequential  addition 
of  solute  gives  a  better  separation  than  batch  loading  for  the  same  number  of  tubes. 

The  partition  coefficients  may  be  stated  at  the  beginning  of  the  simulation  and 
assumed  constant,   or  a  subroutine  JVARY  may  be  written  to  vary  the  partition  co- 
efficients as  a  function  of  solute  concentration.     The  example  of  JVARY  in  this  publi- 
cation shows  the  partition  coefficients  (CK,   XK)  varying  as  a  linear  function  of  total 
amount  of  solute  (T2,   XT2)  in  a  tube.     CK  increases  with  solute  concentration  and 
XK  decreases.     There  is  nothing  fundamental  about  a  linear  variation  and  other  func- 
tions will  be  more  suitable  in  some  systems. 

In  a  countercurrent  distribution,   the  lower  phase  volume  is  not  always  equal 
to  the  cut-off  volume  of  a  tube.     If  the  lower  phase  is  above  the  cut-off  point,   some 
"stationary"  phase  will  be  transported  along  with  the  mobile  phase.     Lower  phase 
can  also  be  held  at  the  air-liquid  interface  and  thus  be  transported.     If  the  lower 
phase  is  below  the  cut-off  point,   some  of  the  mobile  phase  will  be  retained.     Mobile 
phase  can  also  be  retained  on  the  walls  of  the  tubes*.     For  these  reasons  four  volume 
terms  are  needed  to  describe  a  simulation:    total  upper  volume,  total  lower  volume, 
volume  upper  retained,   volume  lower  retained.     These  four  terms  are  input  para- 
meters. 

During  a  simulation,   the  amount  of  overlap  between  two  solute  profiles  may 
be  calculated.     This  option  is  selected  by  turning  on  data  switch  '5'.     The  tubes  in 
the  apparatus  are  separated  into  two  groups  at  the  point  where  the  two  curves  cross. 
The  contents  of  all  tubes  to  the  left  of  the  crossing  point  are  collected  as  the  first 
fraction.     The  contents  of  all  tubes  to  the  right  are  collected  as  the  second  fraction. 
Tubes  with  solute  amounts  less  than  0.  0001  are  ignored.     The  tube  number  and  height 
of  each  peak  are  found.     The  weights  of  both  components  in  each  fraction  are  calcula- 
ted.    The  percent  impurity  of  a  fraction  is  defined  as:    weight  of  lesser  solute/total 
weight  of  both  solutes.     The  total  impurity  is  the  sum  of  the  two  percent  impurities 
and  is  also  in  percent  units.     The  above  computations  are  printed  for  each  transfer. 
However,    if  data  switch  '13'  is  on,   the  impurity  calculations  are  under  the  same  con- 
trol as  the  main  printout. 

Program  I  is  in  two  sections.     The  first  section,   called  SEP62,   reads  the  in- 
put parameters  and  performs  the  simulation,    writing  the  output  profiles  on  the  disk. 
Since  this  calculation  can  be  quite  time-consuming,   the  program  can  be  interrupted 
and  restarted  at  a  later  time.     Data  switch   '3'  will  interrupt  the  calculation  in  this 
manner,    causing  the  restart  information  to  be  saved  on  the  disk.     If  200  tubes  are 
utilized  in  a  simulation,    approximately  350  transfers  an  hour  will  be  performed. 
The  first  10  transfers  are  automatically  printed  and  then  transfers  50,    100,    150,    and 
so  forth  are  printed.     The  operator  can  override  the  automatic  print  control  if  de- 
sired.    While  data  switch  '11'  is  on,    current  transfers  are  printed.     Data  switch  '12' 
on  will  suppress  all  printout.     Switch  '12'  is  especially  useful  when  only  the  output  is 
of  interest. 


*See  reference  in  footnote,    p.    1 


The  output  information  from  SEP62  consists  of  transfer  number,   tube  number, 
and  for  each  of  solute  systems   1  and  2:     total  amount  of  solute,    concentration  in  the 
upper  phase,    concentration  in  the  lower  phase.     Tube  numbering  begins  at  zero;  for 
example,    in  a  200-tube  simulation  the  tubes  will  be  numbered  from  0  to  199.     Trans- 
fer number  also  begins  at  zero,   that  is,    operation  1  is  labeled  transfer  number  0. 

The  second  section  of  Program  I,    called  VEP62,   reads  the  output  profile 
from  the  disk  and  writes  it  on  the  printer.     For  simulations  involving  a  large  number 
of  transfers  only  every  2d,    5th,    and  10th  transfer  may  be  of  interest.     The  program 
reads  one  card  with  the  increment  between  transfer  number.     If  the  card  is  blank, 
every  transfer  will  be  printed.     No  information  will  be  printed  unless  the  total  solute 
eluted  is  greater  than  0.0001.     The  percent  impurity  calculation  is  like  that  in  SEP62. 
This  calculation  is  automatically  printed  unless  data  switch  '5'  is  on. 

The  output  profile  may  be  punched  on  cards  for  plotting  or  other  uses.     To 
punch  the  profile  from  the  first  partition  coefficient  (CK-system),   turn  on  data 
switch  'l1;  from  the  second  partition  coefficient  (XK),    turn  on  data  switch  '2'.     Both 
switches  may  be  on  simultaneously.     The  cards  of  the  XK  output  profile  will  come 
out  in  the  alternate  stacker.     Solute  amounts  less  than  0.0001  will  not  be  punched. 

Card  Input  Required  for  SEP62 


Card  1 


Format 


Number  of 
transfers 


Number  of 
tubes 


215 


Card  2 


Total  upper 
volume 


Total  lower 
volume 


Volume  upper 
retained 


Volume  lower 
retained 


4F10.5 


Card  3 


Kp  first 
solute 

Card  4 


Grams  per 
increment 

first  solute 


Number  of  increments 
first  solute 


2F10.  5,  15 


Grams  per 
K-q  second  increment  Number  of  increments 

solute  second  solute  second  solute 


2F10.5,  15 


Sample  Input  for  SEP62 


500    200 

20.  0 

39.2 

1. 

4 

5.  2 

1.0 

50 

7.  6 

2.0 

35 

38.  6 


Card  Input  for  VEP62 


Card  1 


Increment 
for  printout 


Item 

Partition  coefficient 

Upper  phase  volume 

Lower  phase  volume 

Number  of  tubes 
in  apparatus 

Number  of  increments 
of  solute 

Upper  volume 
transferred 

Lower  volume 
transferred 


Format 


13 


Correlation  of  Systems  of  Symbols 


Program 
CK,   XK 
VUPPR 
VLOW 

NTUBE 

NOFF1,    NOFF2 

VUPUL 

VLPUL 


Reference* 


K 


D 


V 


U 


vn 


V 


UT 


V 


LT 


Summary  of  Switches  Used  in  SEP62 

'0'  continues  with  previously  interrupted  calculation 

'3'  interrupts  calculation 

'5'  calculates  and  prints  percent  impurities 
'11'  prints  every  transfer 
'12'  suppresses  all  printout 
'13'  places  impurities  printout  under  control  of  automatic 

printout  and  subject  to  control  by  switches  '11'  and  '12' 
'14'  causes  batch  rather  than  sequential  loading  of  solute 

Summary  of  Switches  Used  in  VEP62 

•1'  punches  output  profile  of  the  CK  system 
'2'  punches  output  profile  of  the  XK  system 
'5'  suppresses  percent  impurity  calculation 


*See  reference  in  footnote  p.    1. 


PROGRAM  I 

MAINLINE  PROGRAM  SEP62 

♦EXTENDED  PRECISION 

*ONE  WORD  INTEGERS 

*  I OCS ( CARD, 1 132  PR INTER,D I SK, TYPEWRITER) 

AUGUST  5,1969   V.G.  MARTIN 

ARRAY  CONTAINING  CONCENTRATIONS. 

PARTITION  COEFFICIENT   (XK) 

CONCENTRATION  IN  LOWER  PHASE  <XCL01,2) 

CONCENTRATION  IN  UPPER  PHASE  (XCUP1,2) 

ARRAY  FOR  OUTPUT  PROF  I LE . 

TUBE  NUMBER  PLUS  ONE. 

RECORD  NUMBER  IN  FILE  2  -  APPEARS  IN  DEFINE  FILE 

STATEMENT. 

TRANSFER  NUMBER  PLUS  ONE 

LOWER  BOUND  OF  N  IN  TRANSFER  NUMBER  LOOP. 

ARGUMENT  3  IN  SUBROUTINE  MOD.   IT  EQUALS  THE  INTEGER 

PART  OF  ARGUMENT  1  DIVIDED  BY  ARGUMENT  2. 

INDICATOR  FOR  SW I TCH  3. 

1  EXIT  AFTER  PRINTING  OUT  NEXT  COMPLETE  JR-LOOP. 

2  CONTINUE  WITH  CALCULATION. 
UPPER  BOUND  OF  N  IN  TRANSFER  NUMBER  LOOP. 
NUMBER  OF  SOLUTE  INPUTS  (NOFF2) 

SUBSCRIPT  WITHIN  RECORD  OF  FILE  2  =  NEXT  AVAILABLE 
POSITION. 
PRINT  CONTROL 

1  PRINT 

2  DON'T    PRINT 
RECORD    NUMBER    IN    F I LE    1. 
NUMBER    OF    TUBES    IN    APPARATUS. 
TOTAL    SOLUTE    IN    A    TUBE    (XT2) 
GRAMS    OF    SOLUTE    FED    ON    ONE    INPUT    (TINT2) 
LOWER    VOLUME    TOTAL. 
VOLUME    LOWER    PULSED. 
UPPER    VOLUME    TOTAL. 
VOLUME    UPPER    PULSED. 
VOLUME    LOWER    RETAINED. 
VOLUME    UPPER    RETAINED. 

PEAKC    PEAKX 
DIMENSION    FILE2Q04),     BUFFR(IOO) 
DEFINE    Fl LE    1  (9, 320, U, NREC ) 
DEFINE    FILE    2 ( 300, 320, U, KREC ) 
5  FORMATCl    TRANS       TUBE'     13X '  TOTAL '  5X »  CONC .    UP'3X'CONC.     LOW'llX'TOTA 

-L'5X'CONC.    UP'3X'CONC.     LOW1/) 
10  FORMATC    TURN    ON    SWITCH       0       TO    CONTINUE   OLD    RUN'/ 

-'    TURN    ON    SWITCH       5      FOR    PERCENT    IMPURITIES    CALCULATION'/ 
-'    TURN    ON    SWITCH    Ik       FOR    BATCH    LOADING'/) 
15  FORMAT    (215) 

20  FORMAT    UF10.5) 

25  FORMAT    (2F10.5,     15) 

30  FORMAT (I 7,I6,7X,3F12.6,5X,3F12.6) 

35  F0RMAT(//I6,2X, 'TUBES' ) 

kO  F0RMAT(I6,2X'TPANSFERS' ) 


** 

VERSION  OF 

c 

BUFFR  = 

c 

CK 

c 

CL01,2= 

c 

CUP1,2= 

c 

Fl LE2  = 

c 

JR 

c 

KREC   = 

c 

c 

N 

c 

Nl 

c 

N25 

c 

c 

NEXIT  = 

c 

= 

c 

= 

c 

NLIM   = 

c 

NOFF1  = 

c 

NOUT   = 

c 

c 

NPRNT  = 

c 

= 

c 

= 

c 

NREC   = 

c 

NTUBE  = 

c 

T2 

c 

TINT1  = 

c 

VLOW   = 

c 

VLPUL  = 

c 

VUPPR  = 

c 

VUPUL  = 

c 

VLRET  = 

c 

VURET  = 

INTEGER  P 

PROGRAM 


SEP62     (CONT.) 


k5  FORMAT    COUPPER    VO LUME * F 12 . 2  ) 

50  FORMATC     LOWER    V0LUME'F12.2) 

55  FORMATC    UPPER    TRANSFERED'    F8.2) 

60  FORMATC     LOWER    TRANSFERED1    F8.2) 

65  FORMAT    ( /// ' 0  PARTI T I  ON    COEFF .  '  5X ' F EED    WT. * 3X ' NUMBER    OF    TUBES    INITI 

-ALLY'/37X'FI LLED    WITH    SOLUTE1/) 
70  FORMAT(///*0    PARTITION    COEFF .'  5X ' F EED    WT.  '  3X  '  NUMBER    OF    INPUTS1/) 

75  FORMAT(F12.3,F18.3,6X, I  7) 

80  FORMATCOTRANSFER1     15) 

85  FORMATC     PEAK1       LOCATION    ='U'       HEIGHT    = '  F  7 .  k,  7X '  PEAK2       LOCATION    = 

A  '  I  I*  "       HEIGHT    =  '  F  7  .  1+ ) 
90  FORMATC    FRACTION    1 ' 5X ' TUBES ' \k ' - '  I  3, 17X ' FRACT I  ON    2 ' 5X ' TUBES ' I k ' - ' 

Al  3) 
95  FORMATC    WEIGHT    C-COMPONENT    = ' F9 . k, 16X ' WE  I GHT    C-COMPONENT    ='F9.U) 

100  FORMATC  WEIGHT  X-COMPONENT  = ' F 9 . k , 16X 'WE  I GHT  X-COMPONENT  =  '  F 9 . k  ) 
105  FORMATC  PERCENT  IMPURITY  =' F 9 . h , 16X ' PERCENT  IMPURITY  =  '  F 9 . k  ) 
110         FORMAT(21X'TOTAL    IMPURITY    =  '  F  9  .  i+ /  ) 

MAXNO  =  0 

NEXIT  =  2 

WRI TEC  1,10) 

PAUSE  7777 

CALL  DATSW(5,MPURE) 

CALL  DATSW  (l^JSWHO 
C      SWITCH  0  --  OFF-  START  NEW  RUN.   ON--  CONTINUE  OLD  RUN. 

CALL  DATSW(0,MR) 

GO  TO  (115/120)/MR 
115    READ(l'9)N,NLIM,NTUBE,KREC,NOUT,VUPUL,VLPUL,VURET,VLRET,CK,XK,TINT 
ll,TINT2,NOFFl,NOFF2 

NOl  =  NOFF1 

N02  =  NOFF2 

READ(2'KREC)FI  LE2 

KREC  =  KREC  -  1 

Nl  =  N  +  1 

VUPPR  =  VUPUL  +  VURET 

VLOW  =  VLPUL  +  VLRET 

NTUBE  =  NTUBE  -  1 

GO    TO    195 
120         READ(2/15)NLIM/NTUBE 

READ (2, 2  0 ) VUPPR, VLOW, VURET, VLRET 

VUPUL    =    VUPPR-VURET 

VLPUL  =  VLOW  -  VLRET 

READ(2,25)CK,TINT1,N0FF1 

READ(2,2  5)XK,TINT2,NOFF2 

NOl  =  NOFF1 

NO  2  =  NOFF2 

DO  125  KR=1,100 
125    BUFFR(KR)  =  0.0 

DO  130  KR=1,8 
130    WRITE(1'KR)BUFFR 

GO  TO  (135,195),JSW14 
135    CUP1  =  TINT1  /  (VUPPR  +  VLOW/CK) 


PROGRAM  I 


SEP62  (CONT.) 


140 

145 

150 

155 
160 

165 
170 

175 

180 
185 


190 


195 


200 

205 
210 


C 

215 


220 


CLOl  =  CUP1  /  CK 

XCUP1=  T!NT2/(VUPPR  +  VLOW/XK) 


XCLOl  =  XCUP1  / 
IF(NOFFl-NOFF2) 


XK 
140,145,145 


MAXNO 
MINNO 
GO  TO 
MAXNO 
MINNO 
DO  190 


=  NOFF2 
=  NOFF1 
150 

=  NOFF1 
=  NOFF2 
1=1, MAXNO 


IF(I-MINNO)  170,170,155 

IF(NOFFl-NOFF2)  160,165,165 

CUP1    =    0.0 

CLOl    =    0.0 

GO    TO    170 

XCUP1    =    0.0 

XCLOl    =    0.0 

J  =  l 

CALL    MOD    (J,25,N25) 

IF(J)    175,175,180 

J    =    25 

GO    TO    18  5 

N25    =    N25    +    1 

READ(1'N25)BUFFR 

BUFFR  (J)  =  CUP1 

BUFFR  (J+  25)  =  CLOl 

BUFFR  (J+  50)  =  XCUP1 

BUFFR  (J+75)  =  XCLOl 

WRITE(1'N25)BUFFR 

NOFF1  =  0 

NOFF2  =  0 

WRITE(3,35)NTUBE 

WRITE(3,40)NLIM 

WRITE(3,45)VUPPR 

WRITE(3,50)VLOW 

WRITE(3,55)VUPUL 

WRITE(3,60)VLPUL 

GO  TO  (200,205),JSW14 

WRITE (3, 65) 

GO  TO  210 

WRITE(3,70) 

WRITE(3,75)CK,TINT1,N01 

WRITE(3,75)XK,TINT2,N02 

GO  TO  (220, 215), MR 

INITIALIZATION  STEPS  FOR  NEW 

KREC  =  1 

NOUT  =  1 

Nl  =  1 

NTUBE  =  NTUBE  +  1 

ISIGN  =  1 

IF(CK-XK)  225,230,230 


RUN, 


PROGRAM  I 

SEP62  (CONT.) 

225  ISIGN  =  -1 
230  WRITE(3,5) 
C  BEGIN    N    -    LOOP-    ------------------    -BEGIN    N-LOOP 

DO    575    N=N1/NLIM 

Nl  =  0 

N2  =  0 

N3  =  0 

N4  =  0 

KDEL  =  2 

SUMC  =  0.0 

SUMX  =  0.0 

TMAX  =0.0 

XTMAX  =0.0 
C      PRINT  ROUTINE 

GO  TO  (255/235)/NEXIT 
235    IF(N-IO)  260,260,240 

C      SWITCH  12  -  ON--DON'T  PRINT  ANYTHING. 
240    CALL  DATSW  (12,  KR  ) 

GO  TO  (265,245),KR 
C      SWITCH  11  -  ON--  PRINT  EVERYTHING. 
245    CALL  DATSW  (  11,  KR  ) 

GO  TO  (260, 250), KR 
250    KR  =  N 
C      SUBROUTINE  MOD  -  CALCULATES  ARGUMENT  1  MODULO  ARGUMENT  2. 

CALL  MOD  (  KR,  50,  N25  ) 

IF(KR)  255,255,265 
255    WRITE(3,5) 
260    NPRNT  =  1 

GO  TO  2  70 
265    NPRNT  =  2 

C      NOFF  ROUTINE  FOR  TERMINATING  SOLUTE  INPUT. 
270    IF(N-NOFFl)  280,280,275 
275    TINT1  =  0.0 
280    IF(N-NOFF2)  290,290,285 
285    TINT2  =  0.0 
290    JR  =  N  +  MAXNO 

IFCN+MAXNO-NTUBE)  300,300,295 
295    JR  =  NTUBE 
300    KR  =  JR 

CALL  MOD  (KR,25,N25) 

IF(JR-NTUBE)  335,305,305 
C  EFFLUENT  PROFILE  SECTION 
305    KR  =  KR  -  1 

IF(KR)  310,310,315 
310    KR  =  KR  +  25 

NREC  =  N25 

GO  TO  320 
315    NREC  =  N25  +  1 
320    READd'NREOBUFFR 
C      FILE2  CONTAINS 
C  SOLUTE  1  AMT  UPR,  AMT  LOW  ELUTED 


PROGRAM 


SEP62     (CONT.) 


Fl LE2(NOUT) 
Fl LE2CNOUT+2 
Fl LE2CNOUT+5 
Fl LE2  (NOUT+7 
KR  =  NOUT 
CALL  MOD  (KR 
IF(KR)    325,3 

325  WRITEU'KREC 
NOUT    =    0 

330         NOUT    =    NOUT 
GO    TO    1*70 

335  NREC  =  N25  + 
READU'NREC) 

C  RECYCLING    TH 

3U0         KR    =    JR 

CALL   MOD(KR, 
IF(KR-l)    3U5 

31*5         NREC    =    N25    + 
WRITEU'NREC 
READ(1'N25)B 
KR    =    25 
GO    TO    36  5 

350  NREC  =  N25  + 
IF(JR-l)    355 

355         CUP1    =    TINT1 
XCUP1    =    TINT 
CLOi    =    0.0 
XCLOl    =    0.0 
GO    TO    370 

360  WRITEU'NREC 
READU'N25)B 
CUP1  =  BUFFR 
CLOI  =  BUFFR 
XCUP1  =  BUFF 
XCLOl  =  BUFF 
READU'NREC) 
GO    TO    370 

365  CUP1  =  BUFFR 
CLOI  =  BUFFR 
XCUP1  =  BUFF 
XCLOl    =    BUFF 

370  CUP2  =  BUFFR 
CL02  =  BUFFR 
XCUP2  =  BUFF 
XCL02  =  BUFF 
T2  =  CUP1*VU 
XT2  =  XCUP1* 
GO    TO    (375,1* 

375         IFU2-TMAX) 

380         TMAX    =    T2 


SOLUTE    2    AMT    UPR,    AMT    LOW    ELUTED 
=    VUPUL    *    BUFFR(KR) 
6)    =    VLPUL    *    BUFFRUR  +  25) 
2)    =    VUPUL    *    BUFFR(KR+50) 
8)    =    VLPUL    *    BUFFR(KR+75) 

,26,N25) 
25,330 
)FI LE2 


1 
BUFFR 
E  JR  LOOP  BEGINS  HERE. 

25,N25) 
,350,365 

1 
)BUFFR 
UFFR 


,355,360 
/VUPUL 
2 /VUPUL 


)BUFFR 

UFFR 

(25) 

(50) 

R(75) 

RUOO) 

BUFFR 

(KR-1) 

(KR+21*) 

R(KR  +  I*9) 

R(KR+7I*) 

(KR) 

(KR+25) 

R(KR+50) 

R(KR+75) 

PUL+CL01*VLPUL+CUP2*VURET+CL02*VLRET 

VUPUL+XCL01*VLPUL+XCUP2*VURET    +XCL0  2*VLRET 

50),MPURE 

385,380,380 


PROGRAM    I 


SEP62     (CONT.) 

PEAK  C    =    JR-1 
385         IFUT2-XTMAX)    395,390,390 
390         XTMAX    =    XT2 

PEAKX=JR-1 
395    DEL  =  (T2-XT2)*ISIGN 

SUMC  =  SUMC  +  T2 

SUMX  =  SUMX  +  XT2 

GO  TO  (435,400), KDEL 
400    IF(DEL)  1*05, 1*20,1*20 
405    KDEL  =  1 

IF(N4)  410,410,415 
1*10    N2  =  JR-1 

GO  TO  i*  50 
415    N3  =  JR 

N2  =  JR-1 

WTC  =  SUMC  -  T2 

WTX  =  SUMX  -  XT2 

SUMC  =  T2 

SUMX  =  XT2 

GO    TO    1*50 
420         IFU2  +  XT2-.0001)    1*50,450,1*25 
425         IF(N4)    1*30,1*30,450 
1*30         Nl*    =    JR-1 

GO    TO    1*50 
1*35         IF(T2  +  XT2-.0001)    1*40, 1*40,  4  50 
440         I F (Nl)    445,445,450 
445         Nl    =    JR 
450         CALL    JVARY(T2,XT2,CK,XK) 

CL02    =    T2    /    (VLOW+CK*VUPPR) 

CUP2    =    CK    *    CL02 

XCL02    =    XT2    /    (VLOW+XK*VUPPR) 

XCUP2    =    XK    *    XCL02 

BUFFR(KR)    =    CUP2 

BUFFRUR+25)  =  CL02 

BUFFRCKR+50)  =  XCUP2 

BUFFRUR+75)  =  XCL02 

CALL  DATSW  (12, MR) 

GO    TO     (455, 460), MR 
455         NPRNT    =    2 
460         GO    TO    (465, 470), NPRNT 
465         M    =    N    -    1 

KR    =    JR    -    1 

WRITE(3,30)M,KR,T2,CUP2,CLO2,XT2,XCUP2,XCLO2 
470         JR    =    JR    -    1 

I F  (JR)    475,475,340 
475         WRITEd'DBUFFR 

GO    TO    (480,555),MPURE 
480         CALL    DATSW(13,JSW13) 

GO    TO    (485,490),JSW13 
485         GO    TO     (490, 555), NPRNT 
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PROGRAM  I 

SEP62  (CONT.) 

490    M  =  N  -  1 

WRITE(3,80)M 
IF(N4)  495,495,530 
495    IF(ISIGN)  500,520,520 
500    IF(TMAX-.OOOl)  505,505,510 
505    MPURE  =  2 

GO  TO  555 
510    WRITE(3,85)PEAKC,TMAX 

SUM1  =  SUMX 
515    WR!TE(3,90)N1,N2 

WRITE(3,95)SUMC 

WRITE(3,100)SUMX 

FIMP  =  SUM1*100./(SUMC+SUMX) 

WRITE(3,105)FIMP 

GO  TO  550 
520    IF(XTMAX-.OOOl)  505,505,525 
525    WRITE(3,85)PEAKX,XTMAX 

SUM1  =  SUMC 

GO  TO  515 
530    IF(ISIGN)  535,545,51+5 
535    WRITE(3,85)PEAKC,TMAX,PEAKX,XTMAX 

SUM1  =  SUMX 

SUM2  =  WTC 
540    WRITE(3,90)N1,N2,N3,N4 

WRITE(3,95)SUMC,WTC 

WRITE(3,100)SUMX,WTX 

FIMP  =  SUM1  *  100./CSUMX+SUMC) 

GIMP  =  SUM2*100./(WTX+WTC) 

WRITE(3,105)FIMP,GIMP 

FIMP    =    FIMP+GIMP 

GO    TO    550 
54  5         WRITE(3,85)PEAKX,XTMAX,PEAKC,TMAX 

SUM1    =    SUMC 

SUM2    =    WTX 

GO    TO    540 
550         WRITE(3,110)FIMP 
555         GO    TO    (570,560),NEXIT 
C  SWITCH    3    -       INTERRUPT    CALCULATION. 

560         CALL    DATSW    (3,KR) 

GO    TO    (565, 575), KR 
565         NEXIT    =    1 

GO    TO    575 
5  70         WRITE(l,9)N/NLIM/NTUBE/KREC/NOUT/VUPUL,VLPUL/VURET/VLRET,CK/XK/TIN 
1T1,TINT2,N0FF1,N0FF2 

WRITE(2'KREC)FI LE2 

CALL    EXIT 
575         CONTINUE 

N    =    N    -    1 

GO    TO    5  70 

END 
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PROGRAM  I 


SUBROUTI NE  MOD 

♦EXTENDED  PRECISION 
*ONE  WORD  I NTEGERS 
C      V.  G.  MARTIN 

SUBROUTINE  MOD ( NUMB R, INCRE,MANY) 
C      CALCULATE  ARGUMENT  1  MODULO  ARGUMENT  2 

MANY  =  NUMBR  /  I NCRE 

NUMBR  =  NUMBR  -  MANY  *  I NCRE 

RETURN 

END 


SUBROUTINE  JVARY 

*ONE  WORD  INTEGERS 
♦EXTENDED  PRECISION 
C      V.  G.  MARTIN 

SUBROUTINE  J  VARY (T2 , XT2 , CK, XK ) 

C VARY  THE  PARTITION  COEFFICENTS  AS  A  FUNCTION  OF  TOTAL 

C  CONCENTRATION  OF  SOLUTE. 

CK  =  .1077  +  . 0007*  T2 

XK  =  .1077  -  .0007  *  XT2 

RETURN 

END 


MAINLINE  PROGRAM  VEP62 

♦EXTENDED  PRECISION 

♦ONE  WORD  INTEGERS 

♦IOCS (CARD, 11 32  PR  INTER, DISK) 

♦♦    VERSION  OF  MAY  26,  1969     V.G.  MARTIN 

INTEGER  PEAKC,PEAKX 

DIMENSION  Fl LE2(10U) 

DEFINE  F I LE  1(   9  ,  32 0, U , NREC ) 

DEFINE  FILE  2 ( 300, 320 , U, KREC ) 
5      FORMAT  (I 6,2E13.5,2X,2E13.5) 

10     FORMAT ( ' OTRANSF ER' 5X' TOTAL' 7X' UPPER' 10 X1 TOTAL* 7X 'UPPER '//  ) 
15     FORMAT  (13) 
20     FORMAT  (2X,2F13.8) 
25     FORMAT( '0PEAK1   LOCATION  ='IU'   HEIGHT  =  '  F 7  .  k , 7X * PEAK2   LOCATION  = 

-'  U'   HEIGHT  ='F7.U) 
30     FORMATC  FRACTION  1 ' 5X ' TUB ES ' \k ' - ' I  3 , 1 7 X ' F RACT I  ON  2  '  5X ' TUBES ' I k ' - ' 

-I  3) 
35     FORMATC  WEIGHT  C-COMPONENT  = ' F9. h, 16X * WE  I GHT  C-COMPONENT  =  '  F 9 . k  ) 
40     FORMATC  WEIGHT  X-COMPONENT  =  '  F 9 . k , 16X * WE  I GHT  X-COMPONENT  =  '  F 9 . k ) 
45     FORMATC  PERCENT  IMPURITY    ='  F 9 . k, 16X ' PERCENT  IMPURITY    =  '  F 9 . k  ) 
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PROGRAM  I 


VEP62  (CONT.) 

50     FORMAT(21X'TOTAL  IMPURITY  ='F9.4) 
CALL  DATSW  (1,JSW1) 
CALL  DATSW  (2 „ JSW2  ) 
CALL  DATSW(5/MPURE) 
READ(2/15)INCR 
I F ( JSW1  +  JSW2  -  h)    52,  53,  53 

52  READ  (2,15) 

53  IF(INCR)  55,  55,  60 
55     INCR  =  1 

60     WR!TE(3,10) 

READ(l'9)N,NLIM/NTUBE,KREC/NOUT/V/V/V/V/CK,XK 

GO  TO  (80,65),MPURE 
65     ISIGN  =  1 

KDEL  =  2 

IF(CK-XK)  70,75,75 
70     ISIGN  =  -1 
75     SUMC  =  0.0 

SUMX  =  0.0 

TMAX  =0.0 

XTMAX=  0.0 

Nl  =  0 

N2  =  0 

N3  =  0 

Nk    =    0 
80     K  =  NTUBE  -  2 

KREC  =  1 

DO  320  l=NTUBE,N,26 

READ(2'KREC)FI LE2 

DO  315  J=l,26 

K  =  K  +  1 

Kl  =  K  /  INCR 

Kl  =  Kl  *  INCR  -  K 
85  IF(K-N)  90,325,325 
90     TUP  =  Fl  LE2(J) 

TLO  =  Fl LE2CJ+26) 

T2  =  TUP+TLO 

XTUP  =  Fl LE2CJ+52) 

XTLO  =  Fl LE2(J+78) 

XT2  =  XTUP  +  XTLO 

TEST  =  T2  +  XT2 

IF(TEST-.OOOl)  95,95,125 
95     IF(N1)  315,315,100 
100    I F ( N 2  )  105,105,110 
105    N2  =  K  -  1 

GO  TO  315 
110    IF(N3)  315,315,115 
115    I  F  (N£+ )  120,120,315 
120    NU  =  K  -  1 
125    IF(K1)  160,130,160 
130    WRITE(3,5)K,T2,TUP,XT2,XTUP 
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PROGRAM 


VEP62  (CONT.) 

GO  TO  (135, li*5), JSW1 

135    XK  =  K 

IFCT2-1.0E-0U)  11*5,145,140 
140    WRITE(2,20)XK,T2 
145    GO  TO  (150,160),JSW2 
150    XK  =  K 

I F (XT2-1 . OE-OU  )  160,160,155 
155    WRITE(2,20)XK,XT2 

CALL    STACK 
160         GO    TO    (310,165),MPURE 
165         IF(T2-TMAX)    175,170,170 
170         TMAX    =    T2 

PEAKC    =    K 
175         I  FUT2-XTMAX)    185,180,180 
180         XTMAX    =    XT2 

PEAKX  =  K 
185    DEL  =  (T2-XT2)*ISIGN 

SUMC  =  SUMC  +  T2 

SUMX  =  SUMX  +XT2 

GO  TO  (190,230),KDEL 
190    IF(ISIGN)  195,225,225 
195    I F (T2- . 0001 )  200,200,215 
200    I FCN3)  310,310,205 
205    IF  (NO  210,210,310 
210    N4  =  K  -  1 

GO  TO  310 
215    IF(N3)  220,220,310 
220    N3  =  K 

GO  TO  310 
225    IF(XT2-.0001)  200,200,215 
230    IF(DEL)  275,235,235 
235    I F( ISIGN)  270, 240,240 
240    I F (T2-. 0001 )  245,245,260 
245    IF(N1)  310,310,250 
250    IF(N2)  255,255,310 
255    N2  =  K-l 

GO  TO  310 
260    IF(N1)  265,265,310 
265    Nl  =  K 

GO  TO  310 
270    IF(XT2-.0001)  245,245,260 
275    KDEL  =  1 

IF(N2)  280,280,285 
280    N2  =  K  -  1 

N3  =  K 

GO  TO  300 
285    IF(ISIGN)  290,290,305 
290    IFCT2-.0001)  300,300,295 
2  95    N3  =  K 
300    WTC  =  SUMC  -  T2 
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PROGRAM  I 


VEP62  (CONT.) 

WTX  =  SUMX  -  XT2 

SUMC  =  T2 

SUMX  =  XT2 

GO  TO  310 
305    IFCXT2-.0001)  300,300,295 
310    CONTINUE 
315    CONTINUE 
320    CONTINUE 
325    GO  TO  (350,330),MPURE 
330    IF(ISIGN)  340,31*0,335 
335    WRITE(3,25)PEAKC,TMAX,PEAKX,XTMAX 

SUM1  =  WTX 

SUM2  =  SUMC 

GO  TO  34  5 
34  0    WRITE(3,25)PEAKX,XTMAX,PEAKC,TMAX 

SUM1  =  WTC 

SUM2  =  SUMX 
345    WRITE(3,30)N1,N2,N3,N4 

WRITE(3,35)WTC,SUMC 

WRITE(3,40)  WTX, SUMX 

FIMP  =  SUM1  *  100./  (WTC+WTX) 

GIMP  =  SUM2  *  100./  (SUMC+SUMX) 

WRITE(3,45)FIMP,GIMP 

F IMP  =  FIMP  +  GIMP 

WRITE(3,50)FIMP 
350         CALL    EXIT 

END 
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PROGRAM  II 

Program  II  performs  a  counter  cur  rent  simulation  of  one  solute.     Information 
from  the  ternary  diagram  of  this  solute  and  the  two  solvents  used  in  the  distribution 
is  incorporated  in  the  computer  simulation.     Many  of  the  limitations  of  Program  I 
are  thus  overcome.     Program  II  allows  for  solute  volume,   solute -solvent  interaction, 
mutual  solubilities  of  the  solvents,    complex  variation  in  the  partition  coefficient, 
volume  reorganization,    and  isopycnic  phenomena*.     All  data  particular  to  a  certain 
ternary  system  are  isolated  in  subroutines  so  the  mainline  programs  need  not  be  re- 
compiled for  each  separate  system.     Program  II  is  general  for  most  Type  I  and  Type 
II  liquid-liquid  systems,  **  provided  the  equilibria  are  expressed  in  the  following 
manner: 

1.  Triangular  coordinates  are  used  and  the  units  are  expressed  in 
weight  fraction. 

2.  The  solute  is  the  independent  variable. 

3.  The  best  solvent  for  the  solute  is  the  dependent  variable  (Y). 
This  convention  makes  the  tie  line  slopes  AY/ AX    positive. 
Solutropes,   therefore,    cannot  be  processed  with  this  program. 

4.  The  other  solvent  (Z)  is  not  specified  but  can  always  be  cal- 
culated by  Z  =  1  -  X  -  Y,    since  the  weight  fractions  add  up  to 
1.  0. 

5.  Tie  lines  are  described  by  their  slopes  and  Y-intercepts  in  a 
subroutine  SBTIE. 

6.  The  solubility  curves  are  expressed  as  polynomials  in  sub- 
routines CURVU,    CURVL,  KURVU,   KURVL.     The  U  and  L 
symbolize  the  upper  and  lower  curves  on  the  triangular  plot, 
not  necessarily  upper  and  lower  phase. 

The  computation  is  done  by  four  mainline  programs:     TERN1,    TERN2, 
TERN3  perform  the  simulation,    while  PUGET  reads  and  prints  the  output  profile. 
TERN1,    TERN2,    and  TERN3  need  not  be  distinct  segments  on  a  computer  with  larger 
core  capacity.     The  main  computation  is  interruptable  by  switch  '6'. 

If  the  program  is  manually  interrupted  by  switch  '6',  TERN3  is  called  and 
saves  the  data  on  the  disk.     The  next  time  the  program  is  executed,   this  interrupted 
run  will  automatically  continue  until  interrupted  again  with  switch  '6',   terminated 
with  switch  '3',    or  the  desired  number  of  transfers  are  completed.     In  the  latter 
two  instances,    it  is  not  possible  to  restart  an  old  run  and  an  execution  will  result  in 
reading  a  set  of  data  and  starting  a  new  computation. 


'See  reference  in  footnote  p.    1. 

'Treybal,    R.   E.     Liquid  extraction.     1963.     Ed.    2.     McGraw  Hill,    New  York. 
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In  segment  TERN1,   the  input  cards  are  interpreted  and  initial  tube  conditions 
are  established.     Other  necessary  input  parameters  are  also  set.     In  the  case  of  a 
restart,    cards  are  not  read  but  the  information  is  taken  off  the  disk. 

TERN2  performs  the  actual  simulation.     The  total  composition  in  each  tube  is 
calculated.     The  composition  point  is  tested  to  see  if  it  is  in  the  miscible  region 
(block  3).     If  so,   the  two  tie  lines  between  which  the  point  lies  are  found  and  their 
intersection  calculated.     A  new  tie  line  is  drawn  between  this  intersection  and  the 
point  in  question.     The  intersections  of  the  new  tie  line  and  the  solubility  curves  are 
found  (block  4).     These  intersections  are  the  compositions  in  weight  fractions  of  the 
upper  and  lower  phases.     The  densities  of  the  two  phases  are  calculated  and  tested 
for  phase  inversion  (block  5).    The  weight  of  each  phase  is  calculated  by  the  lever 
rule.     The  volumes  are  determined  and  a  transfer  is  made  correcting  for  volumes 
retained  on  the  tube  walls  (block  6). 

PUGET  is  a  mainline  program  for  printing  the  output  profile.     This  program 
also  uses  data  from  the  ternary  diagram.     It  performs  essentially  the  same  set  of 
calculations  as  TERN2. 

Input 

Card  1  Format 

Print  control  KPRNT  array  1615 

transfers  from  KPRNT  (N)   -  KPRNT  (N+l) 
will  be  automatically  printed,    N  =   1,    3,    5,...  15 


Card  2 


Mnemonics  for  the  three  components  3  (A3,    IX) 


Card  3 


Number  of  transfers,   number  of  tubes, 

number  of  inputs  315 


Card  4 


Densities  of:    X,    Y,    Z,   most  dense  solvent, 

least  dense  solvent  5F7.4 


Card  5 


Card  6 


Upper,   lower  volumes  retained  on  the  wall  2F5.  2 

Feed  solution  information  19X,    II,    4F10.4 

NGOTO,   XUP1,    XUP2,   XUP3,    VIN 
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There  are  three  choices  for  specifying 
the  feed  solution 

NGOTO  =  1 

XUP1  is  interpreted  as  grams  of  solute  per  input 

XUP2  is  volume  of  pre -equilibrated  upper  phase  to  be  fed 
NGOTO  =  2 

XUP1  is  weight  of  solute  (X) 

XUP2  is  weight  of  Y 

XUP3  is  weight  of  Z 
NGOTO  -  3 

XUP1  is  weight  fraction  of  X  in  feed  solution 

XUP2  is  weight  fraction  of  Y  in  feed  solution 

XUP3  is  weight  fraction  of  Z  in  feed  solution 

VIN  is  the  volume  of  upper  phase  to  be  fed 

Card  7 

Initial  tube  condition  information  19X,    II,    4F10.4,    215 

NGOTO,  XLOW1,  XLOW2,  VLOW,  VECUT,  NTUBA,  NTUBZ 

There  are  two  choices  for  specifying  initial  tube 
conditions 

NGOTO  =  1 

XLOW1  is  weight  of  solute  (X)  initially  in 
tubes  NTUBA  -  NTUBZ 

VLOW  is  volume  of  lower  phase  in 
tubes  NTUBA  -   NTUBZ 

VECUT  is  cut-off  volume  of  tubes 
NTUBA  -   NTUBZ 

NGOTO  =  2 

XLOW1  is  weight  of  solute  (X)  in  tubes 
NTUBA  -   NTUBZ 


XLOW  is  weight  of  Z  in  tubes  NTUBA    -  NTUBZ 

VLOW  is  weight  of  Y  in  tubes  NTUBA   -  NTUBZ 

VECUT  is  the  cut-off  volumes  of  tubes  NTUBA   -   NTUBZ 

If  option  2  is  selected,   some  tubes  can  be  filled  with  both  phases 
initially. 

Depending  on  the  choice  of  NTUBA  and  NTUBZ,    many  different 
number  7  cards  may  be  read  for  a  single  run.     This  permits  partial 
batch  loading  of  the  instrument,    different  cut-off  volumes  for  every 
tube,   and  other  flexibilities.     However,    small  random  deviations  in 
the  cut-off  volumes  of  tubes  such  as  occur  in  a  real  apparatus  have 
been  shown  to  produce  little  change  in  a  distribution*.     The  com- 
puter will  continue  reading  '7'  cards  until  a  value  of  NTUBZ  cor- 
responding to  the  serial  number  of  the  last  tube  +2  is  encountered. 
A  value  of  2  for  NTUBA  corresponds  to  the  first  tube  in  the  ap- 
paratus.    For  example,    if  tubes  0-199  were  to  be  filled  NTUBA  = 
2,    NTUBZ  =  201. 

Card  8 

08888  15       • 

If  you  wish  to  do  another  run  immediately,    omit  the  number  8 
card  and  insert  new  set  of  data  starting  with  print  control. 

Sample  Input 

0         25         50        50     100      100     150     150    200     200    250     250    275     275     300    300 

CMX  CMY  CMZ 

180  10      100 

.7870  1.4740   .9970  1.4740  .9970 

1.3   0.  0 

1  5.1200        19.4900 

2  0.0  20.0  60.0  39.2  2  4 
1          0.0                                    39.2           39.2                  5          11 


*Eddy,    C.    R.,    Showell,    J.   S.,    and  Martin,    V.   G.     Uses  of  digital  computers 
in  theoretical  analytical  chemistry.     III.     Some  computational  experiments  on  irregu- 
larities in  countercurrent  distribution.     Journal  of  the  American  Oil  Chemists' 
Society  [in  press]. 
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Interpretation  of  Input  Cards 

Card  1  Print  the  following  transfers:     0-25,    50,    100, 

150,    200,    250,    275,    300 

Card  2  Label  the  components  CMX,    CMY,    CMZ 

Card  3  180  transfers,    10  tubes,    100  inputs 

Card  4  Densities:     CMX,    CMY,   CMZ,    most  dense  solvent, 

least  dense  solvent 

Card  5  1.3  milliliters  of  upper  phase  retained  on  the  wall 

Card  6  NGOTO  =  1,    5.  12  grams  of  solute  per  input, 

19.  49  milliliters  of  upper  phase  fed  on  each 
operation 

Card  7-a      Fill  tubes  0-2  with  0.  0  gram  CMX, 

20.  0  grams  CMZ,    60.  0  grams  CMY 
Cut-off  volume  for  tubes  0-2  is  39.  2 

Card  7-b      Fill  tubes  3-9  with  39.  2  milliliters  of  lower 
phase  pre-equilibrated  with  upper  phase 
Cut-off  volume  of  tubes  3-9  is  39.  2 

Card  8  When  simulation  is  done,    stop.     Do  not 

proceed  with  another  run. 

Input  for  PUGET 

The  density  card  (card  4)  is  used  as  input  for  PUGET 

Procedure  for  Preparing  Subroutines 

Experimental  data  for  the  upper  and  lower  solubility  curves  must  be  reduced 
to  (5th  degree  or  lower  order)  polynomial  relationships  between  X  (the  solute)  and  Y 
(the  best  solvent).     In  a  Type  I  system,    where  there  is  only  one  solubility  curve,    an 
artificial  division  into  upper  and  lower  curves  may  be  made  at  the  plait  point.     The 
fit  curves  need  only  be  accurate  descriptions  of  the  system  within  the  concentration 
regions  of  interest  in  the  program.     In  our  acetone-water-chloroform  system, 
(fig.    1),  *  only  the  solid  parts  of  the  solubility  curve  were  fit.     The  area  above  tie 
line  1  was  never  reached  in  the  simulation.     Interpolation  beyond  this  tie  line  would 
be  difficult.     Small  errors  near  the  plait  point  lead  to  large  errors  in  relative  phase 
volumes. 


-Brancker,   A.   V.,    Hunter,    T.    G.  ,    and  Nash,   A.    W.     J.    Phys.    Chem.    44: 
683-698.     1940. 
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Figure  1.  --Ternary  diagram  of  acetone, 
water,  chloroform  system  at 
25°  C. 


Figure  2.  --Ternary  diagram  of  methyl 
palmitate,  hexane,  acetoni- 
trile  system  at  25°  C. 


The  polynomials  of  the  upper  and  lower  curves  are  inserted  in  the  appropri- 
ate subroutines.     The  coefficients  of  the  upper  curve  are  used  in  CURVU  and  KURVU. 
In  CURVU  the  coefficients  are  expressed  as  an  array:     A(l)  Y-intercept,    A(2)  co- 
efficient of  X,   A(3)  coefficient  of  X   ,   and  so  forth.     In  KURVU  the  coefficients  are  in 
standard  form  for  evaluating  a  polynomial  by  Horner's  rule*.     The  coefficients  of  the 
lower  curve  are  similarly  used  in  CURVL  and  KURVL.     CURVU  and  CURVL  in  turn 
call  NWRPH,    which  finds  the  roots  of  a  polynomial  by  the  Newton-Raphson  method*. 
In  CURVU  the  value  of  X2  must  be  chosen  according  to  the  types  of  system.     In  a 
Type  II  (fig.    2)  system  X2  is  set  equal  to  the  X-intercept  of  the  upper  curve**.     In  a 
Type  I  system  X2  is  the  highest  value  of  X  on  the  upper  curve  which  is  of  interest  in 
the  simulation.     In  figure  1  this  corresponds  to  the  point  where  tie  line  1  intersects 
the  upper  curve. 

The  equation  of  each  tie  line  is  calculated.     These  equations  are  used  in  sub- 
routine SBTIE.     TIE(J,  1)  is  the  intercept  of  tie  line  J,    TIE(J,  2)  is  the  slope  of  tie 
line.  J.     The  tie  lines  must  be  ordered  in  the  array  as  shown  on  the  ternary  diagram, 
(fig.    1).     The  first  two  positions  in  the  array  TIE  denote  the  intersections  of  the  solu- 
bility curves  with  the  x-axis.     If  there  are  no  x-intercepts,   the  values  2.  0,    2.  0  are 
inserted  in  the  first  two  positions  of  the  TIE  array.     These  values  are  tested  by  the 
mainline  program  to  determine  whether  the  system  is  Type  I  or  Type  II.     The  last 
two  positions  in  the  TIE  array  denote  the  intersections  of  the  solubility  curves  with 
the  y-axis. 


*McCracken,    D.    D.  ,    and  Dorn,    W.   S.     Numerical  methods  and  FORTRAN 
Programming.     1964.     John  Wiley  &  Sons,    Inc.     New  York. 

**Barford,    R.    A.,    Bertsch,    R.   J.,    and  Rothbart,    H.    L.     J.    Am.    Oil  Chem. 
Soc.     45:     141-143.     1968. 
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In  subroutine  TUBEN,   the  volume  of  upper  phase  to  be  fed  (after  solute  input 
is  discontinued)  must  be  set  at  the  desired  value.     VUP  is  the  name  of  the  variable. 

During  execution  of  TERN2,   the  following  switches  can  be  utilized: 

'3'  ends  program  without  ability  to  restart 

'6'  interrupts  program  with  ability  to  restart 
'11'  every  transfer  printed 
'12'  all  printout  is  suppressed 


During  execution  of  PUGET,   switch  '3'  can  be  used  to  end  printout  when  de- 


sired. 


A  note  on  the  symbols  used  in  the  programs.     These  programs  were  original- 
ly written  for  a  ternary  system  of  methyl  palmitate,   hexane,   acetonitrile.     The  let- 
ters MEP,    HEX,   ACN  or  simply  M,    H,    A  appear  as  part  of  the  variable  names 
throughout  the  programs.     These  letters  in  the  general  case  refer  to  components  X, 
Y,    Z,    respectively. 
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PROGRAM  I  I 

MAINLINE  PROGRAM  TERN1 

*ONE  WORD  INTEGERS 

♦EXTENDED  PRECISION 

*  IOCS (CARD, 11 32  PR  INTER, DISK, TYPEWRITER) 

**  TERN1  -  VERSION  OF  MAY  20,  1968 

C      THIS  PROGRAM  SIMULATES  A  COUNTER  CURRENT  DISTRIBUTION,  MAKING  USE 

C        OF  THE  TERNARY  DIAGRAM  THAT  DESCRIBES  THE  THREE  COMPONENT  SYSTEM 

C        I N  WEIGHT  FRACTIONS. 

C      V.  G.  MARTIN,  R.  A.  BARFORD 

DIMENSION  WTMEP(2  02),WTHEX(2  02),WTACN(202),TI E(20,2) 
DIMENSION  KPRNTU6),  VCUT(202) 
DIMENSION  WDSKK5),  WDSK2(5),  WDSK3(5) 

COMMON  DENSA,DENSH,DENSM,FTACN,FTHEX,FTMEP,GRAMP,DENLT,DENHV 
COMMON  I,  IDENT,  JR,  JROUT,  JRSGL,  JSWCH,  K,  KDISK,  KPRN 
COMMON  KPRNT,  KPRON,  KRITE 

COMMON  L,  N,  NDISK,  NIGHT,  NITE,  NKIM,  NLIM,  NOFF,  NOUT 
COMMON  NREC,  NSTOR,  NTIE,  NTUBA,  NTUBE,  NTUBZ 
COMMON  RHOLO,  RHOUP,  TEM1,  TEM2,  TIE,  TIEB 
COMMON  TIEM,  TRYY,  VCUT,  VECUT,  VI N,  VLACN,  VLHEX,  VLMEP 
COMMON  VLOW,  VLOW1,  VUP,  VUP1,  WALLO,  WALUP,  WDSK1,  WDSK2 
COMMON  WDSK3,  WLMEP,  WTACN,  WTHEX,  WTLOW,  WTMEP,  WTTRA 
COMMON  WTTRH,  WTTRM,  WTTOT,  WTUPR,  WUMEP,  X,  Y,  YTACN,  YTHEX 
DEFINE  FILE  1  ( 9, 320, U, KD I SK ) 
DEFINE  FILE  2  (  200, 53, U, NREC ) 
1  F0RMAT(A3,1X,A3,1X,A3) 

5  FORMAT  ('  COUNTER  CURRENT  DISTRIBUTION  BY  PROGRAM  TERN1,  VERSION 
510F  MAY  20,  196  8'  ) 
10  FORMAT  (1615) 
15  FORMAT  (5F  7.4,215) 

20  FORMAT  Cl   TUBE     WEIGHT  OF       WEIGHT  OF       WEIGHT  TOT     WE 
1IGHT  TOT       WEIGHT  TOT       VOLUME  VOLUME      TRANSFER'/ 

2'   NUMBER     CMX  IN',9X,'CMX  I  N  '  ,  1  2X,  '  CMX  ' ,  1 IX,  '  CMY  '  ,  12X,  '  CMZ,  '  11X 
3, 'UPPER', 1 OX, 'LOWER1, 6 X, ' NUMBER '/ 1 3X, 'U PPER ', 10 X, ' LOWER '/ ) 
35  FORMAT(  19X,  II,  4F10.lt,  215  ) 
kO    FORMAT(//28X, '-INITIAL  CONDITIONS-'/) 
1*1  FORMAT(30X,  'UPPER',  8X,  'LOWER'/) 
4  2  FORMAT (1 OX,  'WE  I GHT ' , A  4 , 5X, F 1 1 . h , F 1 3 . h ) 
1*3  FORMAT(2F5.2) 

kh    FORMATdOX,  'WT  FRCT  ' ,  AU,  4X,  F  11 .  h ,  F  13  .  4  ) 
4  5  FORMATdOX,  '  VO LUME ' ,  9X, Fll . k,  F13.  k ) 

WRITE  (3,5) 
50  CALL  SBTI E  (Tl E,NTI E) 
C      GET  THE  SLOPES  AND  INTERCEPTS  OF  THE  TIE  LINES  FROM  A  SUBROUTINE. 

C BLOCK  1  -  READ  VALUES  FOR  COMPUTATION  FROM  DISK  OR  CARDS. 

100  READ  (I'D  IDENT, DENSM, DENSH, DENSA, NIGHT, NREC, KPRON, KPRN, DENLT, 
ID ENHV, NSTOR, NLIM, NTUBE, NOFF, KRITE, NTUBA, WALUP, WALLO, KPRNT 
IF  (IDENT  -  8)  190,110,110 
C      BEGINNING  OF  NEW  COMPUTATION 
110  READ  (2,10)  KPRNT 
WRITE(3,10)  KPRNT 
IF(KPRNT(l)-8888)  115,210,115 
115  NIGHT  =  0 
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TERN1  (CONT.) 


PROGRAM 


120 
130 


135 


138 


140 


141 


142 


143 


144 


NSTOR  = 
KPRON  = 
KPRN  = 
KRITE  = 
IDENT  = 
EMPTY  F 
DO  120 
WDSKKI 
WDSK2CI 
WDSK3(I 
DO  130 
WRITE  ( 
NREC  = 
READ  BO 
READ(2, 
READ  (2 
WRITE (3 
READ  (2 
DENSM  I 
VARIA 
DENSH  I 
VARIABL 

FOR 

DENSA  I 

THE  0 

WRITE (3 

DENLT  A 

SOL 
READC2, 
WRITE (3 
READ(2, 
GO  TO  ( 
VMEP  = 
WTMEPd 
VIN  =  X 
VUP  =  V 
IF  (DEN 
RHOUP  = 
GO  TO  1 
TEMP=TI 
Tl ECNT1 
Tl ECNTI 
GO  TO  1 
WTUPR  = 
WTHEXQ 
WTACNd 
WTTOT  = 
FUMEP  = 
FUHEX  = 
FUACN  = 
GO  TO  ( 


1 

2 

2 
1 
LE 


I  ) 


2  OF 
1,5 

0.0 
0.0 
0.0 
1,200 
WDSK1, 


PREVIOUS  CONTENTS 


WDSK2,  WDSK3,  NIGHT 


AS  THE  INDEPENDENT 


2 

1 

UNDRY  CONDITIONS  FROM  CARDS. 

1)  CHEM1/CHEM2/CHEM3 

,10)  NLIM,NTUBE,NOFF 

,10)  NLIM,NTUBE,NOFF 

,15)  D ENSM, D ENSH, DENS A, DENHV, DENLT 

S  THE  DENSITY  OF  THE  COMPONENT  EXPRESSED 

BLE  IN  THE  WEIGHT  FRACTION  SYSTEM. 

S  THE  DENSITY  OF  THE  COMPONENT  EXPRESSED  AS  THE  DEPENDENT 

E  IN  THE  WEIGHT  FRACTION  SYSTEM.   IT  IS  THE  BEST  SOLVENT 

DENSM. 
S  THE  DENSITY  OF  THE  COMPONENT  EXPRESSED  AS  A  FUNCTION  OF 
THER  TWO  COMPONENTS  IN  THE  WEIGHT  FRACTION  SYSTEM. 
,15)  DENSM, DENSH, DENSA,DENHV, DENLT 

ND  DENHV  ARE  DENSITIES  OF  LEAST  DENSE  AND  MOST  DENSE 
VENTS  RESPECTIVELY. 
43)  WALUP,WALLO 
,43)  WALUP,WALLO 
35)  NGOTO,XUPl,XUP2,XUP3,VIN 
140,145,150),  NGOTO 
XUP1/DENSM 
)  =  XUP1 
UP2 

IN-VMEP 
SH-DENLT)  4040,  141,  142 

1./ (Tl ECNTI E,2)/DENSH+(1.-TIE(NTIE,2) )/DENSA) 
43 

E(NTIE,2) 
E,2)=TIE(NTI E,l) 
E,l)  =  TEMP 
41 

VUP*RHOUP 
)  =  Tl ECNTI E,2)  *  WTUPR 
)  =  WTUPR-WTHEX(l) 

WTUPR  +  WTMEP(l) 

WTMEP(l)  /  WTTOT 

WTHEX(l)  /  WTTOT 

WTACN(1 )  /  WTTOT 
155,146),  NGOTO 
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PROGRAM  I  I 

TERN1  (CONT.) 

14  5  WTMEP(l)  =  XUP1 

WTHEX(l)  =  XUP2 

WTACN(l)  =  XUP3 

WTTOT  =  WTMEP(l)  +  WTHEX(l)  +  WTACN(l) 

GO  TO  144 
146  RHOUP  =  1.  /  (FUMEP/DENSM  +  FUHEX/DENSH  +  FUACN/DENSA  ) 

VI N  =  WTTOT  /  RHOUP 

GO  TO  155 
150  FUMEP  =  XUP1 

FUHEX  =  XUP2 

FUACN  =  XUP3 

RHOUP  =  1./ (FUMEP/DENSM  +  FUHEX/DENSH  +  FUACN/DENSA) 

WTTOT  =  VI N  *  RHOUP 

WTMEP(l)  =  FUMEP  *  WTTOT 

WTHEX(l)  =  FUHEX  *  WTTOT 

WTACN(l)  =  FUACN*WTTOT 
C      FILLING  THE  TWO  HUNDRED  TUBES  AND  SPECIFYING  THEIR  CUT  OFF  VOLUMES 
C        IS  ACCOMPLISHED  BY  STATEMENTS  140  TO  150. 

155  READ(2,35)  NGOTO , X LO 1 , XL02 , VLOW, VECUT, NTUBA,  NTUBZ 
GO  TO  (156, 157), NGOTO 

156  GRAMP  =  XLOl 

VLOW    =    VLOW-GRAMP/DENSM 

RHOLO    =    l./(TIE(NTIE,l)    /DENSH+ ( 1 . -T I E ( NT  I E, 1 ) ) / DENSA ) 

WTLOW    =    VLOW*RHOLO 

WTTOT    =    GRAMP    +    WTLOW 

YTHEX    =    TIE(NTIE,1)    *    WTLOW 

YTACN    =    WTLOW-YTHEX 

FLMEP    =    GRAMP/WTTOT 

FLHEX    =    YTHEX/WTTOT 

FLACN    =    YTACN/WTTOT 

GO    TO    158 

157  GRAMP  =  XLOl 
YTHEX  =  VLOW 
YTACN=    XL02 

WTTOT=GRAMP+YTHEX+ YTACN 
FLACN=YTACN/WTTOT 
FLHEX=YTHEX/WTTOT 
FLMEP=GRAMP/WTTOT 

158  DO    160     I     =    NTUBA, NTUBZ 
WTMEP( I )    =    GRAMP 
WTHEXd  )    =    YTHEX 
WTACNd  )    =    YTACN 

160  VCUT( I )    =    VECUT 

I F(NTUBZ-NTUBE-l)    155,     161,    4040 

161  I  F(DENSH-DENLT)    4040,     170,     165 
165    TEMP=TI ECNTI  E,l) 

Tl  E(NTIE,1)    =    TIE(NTI E,2) 
TIE(NTIE,2)=TEMP 
170    WRITE(3,40) 
WRITE(3,41) 
WRITE(3,42)    CHEM1,WTMEP(1), GRAMP 
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TERN1     (CONT.) 


PROGRAM    I  I 


W R I TE ( 3 , h 2 )    CH EM2 , WTHEX ( 1 ) , YTH EX 

WRITE (3, 42)    CHEM3,WTACN(1),YTACN 

WRITE(3,44)CHEM1,FUMEP,FLMEP 

WRITE(3,i|4)    CHEM2/FUHEX/FLHEX 

WRITE  (3,  1*4)    CHEM3,FUACN,FLACN 

WRITE  (3,  45)    V  INFLOW 

GO    TO    2  00 

CONTINOE  PREVIOUS  CO 


190 


200 

210 
i+  QUO 


READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

NREC 

WRITE 

CALL 

CALL 

STOP 

END 


(1'2 
(1'3 
(l'lt 

(l'5 
(1'6 
(1'7 
(1'8 


4PUTATION. 
=1,101) 
=102,202) 
=1,101) 
=102,202) 
=1,101) 
=102,202) 
=  1,101) 
=  102,202) 


GET  VARIABLES  FROM  THE  DISK 


)  (WTMEPd  ), 

)  (WTMEPd  ), 

)  (WTHEX(I), 

)  (WTHEXd  ), 

)  (WTACNC I ), 

)  (WTACNC I ), 

)  (VCUT(I),I 

(1'9)  (VCUTd  ),     I     =    102,202) 
(2 INREC)WDSK1,WDSK2,WDSK3,NDISK 

=    NREC-1 
(3,     20) 
LINK    (TERN2) 
EXIT 

l+OUO 


MAI NL I NE  PROGRAM  TERN2 


*ONE  WORD  I  NT 
•EXTENDED  PRE 
*  I OCS ( 1132  PR  I 
**TERN2  VERSI 
C      V.  G.  M 
REAL  LI 
D I  MENS  I 
DIMENSI 
D I  MENS  I 
DIMENSI 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
DEF INE 
2  0  FORMAT( 

C BLOCK  2 

2  00  NKIM  = 
DO  880 


EGERS 

CISION 

NTER,DISK) 

ON  OF  MAY  20,  1968 

ARTIN,  R.  A.  BARFORD 

NEA,  LINEH 

ON  WTMEP(202),WTHEX(2  02),WTACN(2  02),TI E(20,2) 

ON  TEMPO(3) 

ON  KPRNT(16),  VCUT(202) 

ON  WDSKK5),  WDSK2(5),  WDSK3(5) 

DENSA,DENSH,DENSM,FTACN,FTHEX,FTMEP,GRAMP,DENLT,DENHV 

I,  IDENT,  JR,  JROUT,  JRSGL,  JSWCH,  K,  KDISK,  KPRN 

KPRNT,  KPRON,  KRITE 

L,  N,  NDISK,  NIGHT,  NITE,  NKIM,  NLIM,  NOFF,  NOUT 

NREC,  NSTUR,  NTIE,  NTUBA,  NTUBE,  NTUBZ 

RHOLO,  RHOUP,  TEM1,  TEM2,  TIE,  TIEB 

TIEM,  TRYY,  VCUT,  VECUT,  VIN,  VLACN,  VLHEX,  VLMEP 

VLOW,  VLOW1,  VUP,  VUP1,  WALLO,  WALUP,  WDSK1,  WDSK2 

WDSK3,  WLMEP,  WTACN,  WTHEX,  WTLOW,  WTMEP,  WTTRA 

WTTRH,  WTTRM,  WTTOT,  WTUPR,  WUMEP, 

Fl LE  2(200, 53, U, NREC) 

1X,I6,7E15.6,I6) 

-    BEGIN    TRANSFER    NUMBER    LOOP. 
N LI  i-1+l 
N=NSTOR,NKIM 
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X,     Y,    YTACN,     YTHEX 


BEGIN    N    LOOP, 


PROGRAM  I  I 

TERN2  (CONT.) 

C      AUTOMATIC  PRINT  CONTROL. 
GO  TO  (210,230), KPRON 
210  I F(N-KPRNT(KPRN)-2)  300,300,220 
220  KRITE  =  2 

KPRN  =  KPRN  +  1 
KPRON  =  2 
GO  TO  300 
230  IF  (KPRN  -  17)  240 , 300 , 4040 
240  IF(N-KPRNT(KPRN)-2)  300,250,250 
250  KRITE  =  1 

KPRN  =  KPRN  +  1 
KPRON  =  1 

C BLOCK  3  -  BEGIN  TUBE  NUMBER  LOOP. 

300  JR  =  N+NTUBA-2 

IF  (JR  -  NTUBE  -  1)  320,320,310 
310  JR  =  NTUBE  +  1 
320  IF  (JR-1)  4040,  650,  325 
C      BEGIN  MAIN  COMPUTATION. 

325  IF  (WTMEP(JR)-(.lE-08)  )  326,  327,  327 
32  6  WTMEP(JR)  =0.0 

327  WTTOT  =  WTMEP(JR)  +  WTH EX ( J R ) +WTACN ( JR ) 
FTMEP  =  WTMEP(JR)  /  WTTOT 
FTHEX  =  WTHEX(JR)  /  WTTOT 
FTACN  =  WTACN(JR)  /  WTTOT 
TEMPO(l)  =  WTMEP(JR) 
TEMPO(2)  =  WTHEX(JR) 
TEMPO(3)  =  WTACN(JR) 
C      TEST  TO  SEE  IF  THE  COMPOSITION  IS  ON  ONE  OF  THE  AXES,  I.E.  BINARY 
IF(FTMEP)  4040,330,360 
330  IF(FTHEX-TIE(NTIE,1))520,520,340 
340  I F (FTHEX-TI E(NT1 E, 2 ) )  350,520,520 
C      TIE(l,l),TIE(l,2),TIE(NTiE,l),  AND  TIE(NTIE,2)  DO  NOT  CONTAIN  THE 
C        SLOPE  AND  INTERCEPT  AS  IN  THE  REST  OF  THE  TIE  LINE  MATRIX,  BUT 
C        CONTAIN  THE  INTERSECTIONS  OF  THE  SOLUBILITY  CURVES  WITH  THE 
C  X  AND  Y  AXES  RESPECTIVELY.   IF  NO  X-INTERCEPTS  EXIST,  ENTER 

C  2.0  FOR  TIE(1,1)  AND  TIE(1,2) 

350  FUMEP  =  0.0 
FLMEP  =  0.0 
FUHEX  =  TIE(NTIE,2) 
FLHEX  =  Tl E(NTI E,l) 
GO  TO  500 
370  IF(FTMEP-TIE(1,D)  520,520,375 

360  IF(FTHEX)  4040,  361,  385 

361  IF(TIE(1,1)-1.0)  370,  370,  520 
3  75  IF(FTMEP-TIE(1,2))380,52  0,52  0 
380  FUHEX  =  0.0 

FLHEX  =  0.0 

FUMEP  =  Tl E(l,2) 

FLMEP  =  TIE(1,1) 

GO  TO  500 
C      TEST  TO  SEE  IF  THE  COMPOSITION  OF  THE  TUBE  IS  IN  THE  MISCIBLE 
C        REGION. 
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PROGRAM    I  I 

TERN2     (CONT.) 

385    CALL    KURVL    (FTMEP,Y) 
C  ANY    POINT    BELOW    Y,    WHICH    IS    ON    THE    LOWER    SOLUBILITY    CURVE,     IS    IN 

C  THE    M I  SCI BLE    REGION. 

IF(Y-FTHEX)    390,     520,520 
390    CALL    KURVU(FTMEP,Y) 
C  ANY    POINT    ABOVE    Y,    WHICH    IS    ON    THE    UPPER    SOLUBILITY    CURVE,     IS    IN 

C  THE    Ml  SCI BLE    REGION. 

IF(FTHEX-Y)     395,520,520 
395    J    =    1 

C BLOCK    4    -    CALCULATE    THE    SLOPE    AND    INTERCEPT    OF    THE    TIE    LINE    FOR 

C  TUBE    JR. 

U00    J=J+1 

TRYY  =  TIE(J,1)  +  TIE(J,2)*FTMEP 
TEMP  =  FTHEX  -  TRYY 
I F(ABS(TEMP)-.001)  460,405,405 
405  IF  (FTHEX  -  TRYY)  4 10, 460 , 4 70 
C      BRANCH  TO  STATEMENT  410  WHEN  THE  FIRST  TIE  LINE  IS  REACHED  BELOW 
C        WHICH  THE  COMPOSITION  LIES. 
410  IF  (J-2)  4040,  415,  430 
415  IF  (TIE(1,1)-1.0)  420,  4050,  4050 
C      STATEMENT  4040  IS  A  CATCHALL  FOR  EVERY  SEEMINGLY  IMPOSSIBLE 
C        SITUATION. 

420  X  =  -Tl E(J,1)/TIE(J,2) 
GO  TO  4  40 
C      CALCULATE  THE  INTERSECTION  OF  TIE  LINES  J  AND  J-l. 
430  L  =  J-l 

TEMP=TIE(L,2)-TIE(J,2) 
X  =  (TIE(J,1)  -  TIE(L,1))/TEMP 
440  Y  =  Tl E(J,1)  +  Tl E(J,2)*X 
C      CALCULATE  A  NEW  Tl E  LINE. 
TEMP=FTMEP-X 
TIEM=(FTHEX-Y)/TEMP 
TEMP=TI EM*FTMEP 
Tl EB=FTHEX-TEMP 
C      FIND  THE  INTERSECTION  OF  THE  NEW  TIE  LINE  WITH  THE  UPPER  AND  LOWER 
C        SOLUBILITY  CURVES.  TO  GET  THE  WEIGHT  FRACTIONS  IN  THE  UPPER  AND 
C        LOWER  PHASES. 

450  CALL  CURVU(TI EM,TI EB, FTMEP, FUMEP, FUHEX  ) 
CALL  CURVL(TIEM,TIEB,FTMEP,FLMEP,FLHEX) 
C      CURVL  AND  CURVU  WILL  IN  TURN  CALL  SUBROUTINE  NWRPH  IF  THE 
C        SOLUBILITY  CURVES  ARE  EXPRESSED  BY  POLYNOMIALS  HIGHER  THAN  THE 
C        SECOND  DEGREE.  NWRPH  FINDS  THE  INTERSECTION  OF  TWO  FUNCTIONS  BY 
C        THE  TRIAL  AND  ERROR  METHOD. THIS  IS  NOT  AN  EXACT  SOLUTION  AND  FOR 
C        THIS  PROGRAM  THE  FIRST  VALUE  FOUND  WITHIN  .0001  OF  THE  ACTUAL 
C        INTERSECTION  IS  RETURNED  AS  THE  CORRECT  VALUE. 
GO  TO  500 
4  60  TIEM  =  Tl E(J,2) 
Tl EB  =  Tl E(J,1) 
GO  TO  4  50 
470  I F ( J-NTI E  +  l )  400,480,4040 
C       BRANCH  TO  480  OCCURS  WHEN  THE  POINT  LIES  BETWEEN  THE  UPPER-MOST 
C        TIE  LINE  AND  THE  Y-AXIS. 
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PROGRAM    I  I 

TERN2    (CONT.) 

U80    X    =    0.0 

GO    TO    kkO 

C BLOCK    5    -    CALCULATE    THE    WEIGHT    OF    THE    UPPER    AND    LOWER    PHASE. 

500    FUACN=1.-FUMEP-FUHEX 

FLACN=1.-FLMEP-FLHEX 
C  CALCULATE    DENSITIES    OF    UPPER    AND    LOWER    PHASES. 

TEMP=FUMEP/DENSM+FUHEX/DENSH+FUACN/DENSA 

RHOUP=l./TEMP 

TEMP=FLMEP/DENSM+FLHEX/6ENSH+FLACN/DENSA 

RHOLO=l./TEMP 

IF(RHOLO    -    RHOUP)    510,520,530 
C  BRANCH    TO    510    OCCURS    IF    CONCENTRATIONS    ARE    ABOVE    THE    ISOPYCNIC 

C  REGION. 

510    TEMP    =    FLMEP 

FLMEP    =    FUMEP 

FUMEP    =    TEMP 

TEMP    =    FLHEX 

FLHEX    =    FUHEX 

FUHEX    =    TEMP 

TEMP    =    FLACN 

FLACN    =    FUACN 

FUACN    =    TEMP 

TEMP    =    RHOUP 

RHOUP    =    RHOLO 

RHOLO    =    TEMP 

GO    TO    530 

OR    IF     IN    A 


c 

BRANCH    TO    520    OCCURS    IF     IN    THE    ISOPYCNIC    REGION, 

c 

MISCIBLE    REGION. 

520 

FLMEP    =    FTMEP 

FUMEP    =    FLMEP 

FLHEX    =    FTHEX 

FUHEX    =    FLHEX 

FLACN    =    FTACN 

FUACN    =    FLACN 

TEMP=FLMEP/DENSM+FLHEX/DENSH+FLACN/DENSA 

RHOLO=l./TEMP 

WTLOW    =    WTTOT 

GO    TO    6  00 
530    AC    =    FTMEP    -    FLMEP 

BC    =    FTHEX    -    FLHEX 

LINEA    =    SQRT(AC*AC    +    BC*BC) 

AC    =    FUMEP    -    FTMEP 

BC    =    FUHEX    -    FTHEX 

LINEH    =    SQRT(AC*AC    +    BC*BC) 

TEMP    =    LINEA    +    LINEH 

LINEH    =    LINEH    /    TEMP 

WTLOW    =    LINEH    *    WTTOT 

C BLOCK    6    -    CALCULATE    THE    WEIGHT    OF    EACH    COMPONENT    THAT    IS    TO 

C  TRANSFERED    TC    TUBE    JR+1,    AND    WHAT    REMAINS    IN    TUBE    JR. 

600    WTUPR    =    WTTOT    -    WTLOW 

VLOW    =    WTLOW/ RHOLO 

VUP      =    WTU PR/ RHOUP 
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TERN2  (CONT.) 

IF(VUP-WALUP)  1601,  1601,  602 
1601  TEMP  =  WALUP 
WALUP  =  VUP 
VUP1  =  VUP  -  WALUP 
WALUP  =  TEMP 
GO  TO  1604 

602  VUP1  =  VUP  -  WALUP 
1604  VLOW1  =  VLOW  -  WALLO 

IF  (  VLOW1  +  VUP1  -  VCUT(JR))  603,  603,  605 

603  TEM1  =  0.0 
GO  TO  6  35 

605    IF    (    VCUT(JR)    -    VLOW1)    610,    620,    630 
C  THIS    TESTS    WHETHER    THE    INTERFACE    IS    ABOVE    OR    BELOW    THE    SIDE    ARM, 

610    TEM1    =    VUP1/VUP*WTUPR 

TEM2    =    (VLOWl-VCUT(JR)  )/VLOW*WTLOW 

GO    TO    6  40 
620    TEM1    =    VUP1/VUP*WTUPR 

TEM2    =    0.0 

GO    TO    640 
630    TEM1    =    (VUP1-VCUT(JR)+VL0W1)/VUP*WTUPR 
6  35    TEM2    =0.0 
640    WTTRM    =    FUMEP*TEM1    +    FLMEP*TEM2 

WTTRH    =    FUHEX*TEM1    +    FLHEX*TEM2 

WTTRA    =    FUACN*TEM1    +    FLACN*TEM2 

JRSG  L    =    2 

WTMEP(JR)    =    WTMEP(JR)    -    WTTRM 

WTHEX(JR)    =    WTHEX(JR)-WTTRH 

WTACN(JR)    =    WTACN(JR)-WTTRA 

GO    TO    660 
650    WTTRM    =   WTMEP(l) 

WTTRH    =    WTHEX(l) 

WTTRA    =    WTACN(l) 

JRSGL    =    1 
660    WTMEPCJR+1)    =    WTMEP(JR+1)    +    WTTRM 

WTHEXCJR+1)    =    WTHEX(JR+1)    +    WTTRH 

WTACN(JR+1)    =    WTACN(JR+1)    +    WTTRA 

C BLOCK    7    -    PRINT    CONTROL 

70  5    CALL    DATSW(12,JSWCH) 

GO  TO  (74  5,710),JSWCH 
710    CALL    DATSW    (11,JSWCH) 

GO  TO  (730,720), JSWCH 
720  GO  TO  (730,745),KRITE 
730  GO  TO  (800, 740), JRSGL 
740    WUMEP    =    FUMEP*WTUPR 

WLMEP    =    FLMEP*WTLOW 

NOUT    =    N-2 

JROUT    =    JR-2 

WRITE (3,  20)  J  ROUT, WUMEP, WLMEP, TEMPO, VUP, VLOW, NOUT 

745  I F(JR-NTUBE)800,800, 750 

750  I F(WTTRM+WTTRH+WTTRA)  4040,800,760 

760  IF(NIGHT)  4040,770,780 

770  NDISK  =  N 
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PROGRAM 


TERN2  (CONT.) 

780  NIGHT  =  NIGHT  +  1 
WDSK1(NIGHT)=WTTRM 

WDSK2CNIGHT)  =  WTTRH 

WDSK3  (NIGHT)=WTTRA 

I F  (NIGHT  -  5)  800, 790,4040 
790  WRITE(2'NREC)WDSK1,WDSK2,WDSK3,NDISK 

NIGHT  =  0 

C-- BLOCK  8  -  RECYCLING  OF  THE  JR-LOOP. 

800  CALL  DATSW(3, JSWCH)  ' 

GO  TO  (810, 820), JSWCH 
810  IDENT  =  8 

GO  TO  900 
820  IF(JR-l)  4040,840,830 
830  JR  =  JR-1 

GO  TO  32  0 
C      END  OF  JR  LOOP.  END  OF  JR-LOOP, 

840  IF(N-NOFF)  860,850,860 
850  CALL  TUBEN 
C      SUBROUTINE  TUBEN  CONTROLS  SOLUTE  INPUT  TERMINATION. 
860  CALL  DATSW(6,JSWCH) 

GO  TO  (870,880), JSWCH 
870  NSTOR  =  N+l 

GO  TO  9  00 
880  CONTINUE 
C      END  OF  N-LOOP.  END  OF  N-LOOP. 

IDENT  =  9 
900  CALL  LINK  (TERN3) 
4040  STOP  4040 

4050  WRITE(3,4051) 

4051  FORMAT( 'OCONC.  IS  BEYOND  LAST  TIE  LINE.') 

4052  GO  TO  4040 
END 

MAINLINE  PROGRAM  TERN3 

*ONE  WORD  INTEGERS 

♦EXTENDED  PRECISION 

*  IOCS (CARD, 11 32  PR  INTER, DISK) 

**  TERN3  -  VERSION  OF  MAY  20,  1968 

C      V.  G.  MARTIN,  R.  A.  BARFORD 

DIMENSION  WTMEP(2  02),WTHEX(2  02),WTACN(2  02),TI E(20,2) 

DIMENSION  KPRNTQ6),  VCUT(202) 

DIMENSION  WDSKK5),  WDSK2C5),  WDSK3(5) 

COMMON  DENSA,DENSH,DENSM,FTACN,FTHEX,FTMEP,GRAMP,DENLT,DENHV 

COMMON  I,  IDENT,  JR,  JROUT,  JRSGL,  JSWCH,  K,  KDISK,  KPRN 

COMMON  KPRNT,  KPRON,  KRITE 

COMMON  L,  N,  NDISK,  NIGHT,  NITE,  NKIM,  NLIM,  NOFF,  NOUT 

COMMON  NREC,  NSTOR,  NTIE,  NTUBA,  NTUBE,  NTUBZ 

COMMON  RHOLO,  RHOUP,  TEM1,  TEM2,  TIE,  TIEB 

COMMON  TIEM,  TRYY,  VCUT,  VECUT,  VI N,  VLACN,  VLHEX,  VLMEP 
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TERN3  (CONT.) 

COMMON 
COMMON 
COMMON 
DEFINE 
DEFINE 
25  FORMATC 
C BLOCK  9 

900  IFCNIGH 

910    NITE    = 
DO    920 
WDSKKJ 
WDSK2(J 

920  WDSK3CJ 
WRITE (2 
NREC  = 

930  WRITEd 
1DENHV,N 
WRITE (3 
WRITECl 
WRITEd 
WRITECl 
WRITECl 
WRITEd 
WRITEd 
WRITECl 
WRITEd 
IFCIDEN 
50  CALL  LI 

9^0  CALL  EX 
i+OUO  STOP  40 
END 


VLOW, 
WDSK3, 
WTTRH, 
Fl  LE  1 
Fl  LE  2 
//17I7 
-  SAV 
T)I+0U0 
NIGHT+ 
J  =  Nl 
)  =  0. 
)=  0.0 
)  =  0. 
'NREC) 
NREC  - 
•l)  ID 
STOR,N 
,25)ID 
(W 
(W 
(W 

cw 

(W 
(W 

cv 

(V 
)  9 
(TE 


VLOW1,  VUP,  VUP1,  WALLO,  WALUP,  WDSK1,  WDSK2 


WLMEP,  WTACN,  WTHEX, 
WTTRM,  WTTOT,  WTUPR, 

(9/320/U/KDISK) 

(200/53/U/NREC) 

) 

E  DATA  ON  DISK. 

,930,910 

1 

TE,5 

0 


0 


WTLOW,  WTMEP,  WTTRA 
WUMEP,  X,  Y,  YTACN,  YTHEX 


•2) 

'3) 

'4) 

'5) 

•6) 

•7) 

•8) 

'9) 

T-8 

NK 

IT 

40 


WDSK1,WDSK2,WDSK3,NDISK 

1 
ENT,DENSM,DENSH,DENSA,NIGHT,NREC,KPRON,KPRN,DENLT, 
LIM,NTUBE,NOFF,KRITE,NTUBA  ,  WALU  P,  WALLO,  KPRNT 
ENT, KPRN,KPRON, NIGHT, NREC,  KR  I TE, NSTOR, NL I M, NOFF, NTUBE 


TMEPC 

TMEP( 

THEX( 

THEX( 

TACN( 

TACN( 

CUT( 

CUT( 

40,  50 

RN1) 


I  ) 
I  ) 


, 1=1,101) 
,  1=102,202) 
,  1=1,101) 
, 1=102,202) 
,1=1,101) 
, 1=102,202) 
I  =  1,101) 

I  =  102,202) 
940 


SUBROUTINE  TUBEN 

•ONE  WORD  I NTEGERS 
•EXTENDED  PRECISION 
**  TUBEN  -  VERSION  OF 


OCTOBER  27,  1967 


V.  G.  MARTIN 

SUBROUTINE  TUBEN 

DIMENSION  WTMEP(202),WTHEX(2  02),WTACN(2  02),TI E(20,2) 

DIMENSION  KPRNT(16),  VCUT(202) 

DIMENSION  WDSK1[5),  WDSK2(5),  WDSK3(5) 

COMMON  DENSA,DENSH,DENSM,FTACN,FTHEX,FTMEP,GRAMP,DENLT,DENHV 

COMMON  I,  IDENT,  JR,  JROUT,  JRSGL,  JSWCH,  K,  KDISK,  KPRN 

COMMON  KPRNT,  KPRON,  KRITE 

COMMON  L,  N,  NDISK,  NIGHT,  NITE,  NKIM,  NLIM,  NOFF,  NOUT 

COMMON  NREC,  NSTOR,  NTIE,  NTUBA,  NTUBE,  NTUBZ 

COMMON  RHOLO,  RHOUP,  TEM1,  TEM2,  TIE,  TIEB 

COMMON  TIEM,  TRYY,  VCUT,  VECUT,  VI N,  VLACN,  VLHEX,  VLMEP 

COMMON  VLOW,  VLOW1,  VUP,  VUP1,  WALLO,  WALUP,  WDSK1,  WDSK2 
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TUBEN    (CONT.) 

COMMON    WDSK3,    WLMEP,    WTACN/    WTHEX,    WTLOW,    WTMEP,    WTTRA 

COMMON    WTTRH,    WTTRM,    WTTOT,    WTUPR,    WUMEP,    X,    Y,    YTACN,    YTHEX 
C  VUP    IS    THE    VOLUME    OF    PRE-EQU I L I BRATED    UPPER    PHASE    TO    BE    FED 

C  AFTER    SOLUTE    INPUT    IS    SHUT    OFF. 

VUP=20.35 

IF(DENSH-DENLT)  1*040,  101,  102 

101  TEMP=TIE(NTI E, 2 )/DENSH+ ( 1 . -Tl E ( NT  I E,2))/DENSA 
GO  TO  10  3 

102  TEMP=TIE(NTIE/2) 
TIE(NTIE,2)=TIE(NTIE,1) 
Tl E(NTI E/1)=TEMP 

GO  TO  101 

103  RHOUP=l./TEMP 
WTUPR    =    VUP*RHOUP 
WTMEP(l)    =0.0 

WTHEX(l)    =    TIE(NTI E/2)*WTUPR 
WTACN(l)    =    WTUPR    -    WTHEX(l) 
IF(DENSH-DENLT)    1*040,     106,     105 

105  TEMP=TI E(NTI E,2) 
TIE(NTIE/2)=TIE(NTIE/1) 
Tl E(NTIE/1)=TEMP 

106  CONTINUE 
RETURN 

1*040    STOP    4040 
END 


SUBROUTINE  NWRPH 

*ONE  WORD  INTEGERS 

♦EXTENDED  PRECISION 

**   NWRPH  VERSION  OF  SEPTEMBER  21,  1967      V.G.MARTIN 

SUBROUTINE  NWRPH  ( X , CO EFF, X 1, X2  ) 

DIMENSION  COEFF(6) 

AO  =  COEFF(l) 

Al  =  COEFF(2) 

A2  =  COEFF(3) 

A3  =  COEFFU) 

Ak    =  COEFF(5) 

A5  =  COEFF(6) 
1X3=  (XI  +  X2)  /  2. 

X  =  X3 

Y  =((((A5*X  +  AI*)*X  +  A3)*X  +  A2)*X  +  A1)*X  +  A0 

IF(ABS(Y)  -  .0001)  5,2,2 

2  IF  (Y)  3,5,1* 

3  X2  =  X3 
GO  TO  1 

l*  XI  =  X3 

GO  TO  1 
5  RETURN 

END 
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PROGRAM  I  I 


THE  NEXT  FIVE  SUBROUTINES  DESCRIBE  THE  TYPE  I  SYSTEM  IN  FIG.(l) 


SUBROUTINE  CURVU   TYPE  I 

♦EXTENDED  PRECISION 
*ONE  WORD  I NTEGERS 
C      R.  A.  BARFORD 

SUBROUTINE  CURVU ( XM, XB, FTX, X, Y ) 

DIMENSION  A(6) 

A(l)  =  0.99itUU885-XB 

A(2)=-1.16U6217-XM 

A(3)=0. 91691*23 

AU)  — 1. 1*991*1+92 

A(5)=0.0 

A(6)=0.0 

X1=FTX 

X2=0.570 

CALL    NWRPH(X,A,X1,X2) 

Y=XM*X+XB 

RETURN 

END 


SUBROUTINE    CURVL      TYPE    I 


*ONE    WORD 
♦EXTENDED 


INTEGERS 
PRECISION 


R.    A.    BARFORD 

SUBROUTINE  CURVL( XM, XB , FTX, X, Y ) 

DIMENSION  A(6) 

005733911+-XB 

09026U683-XM 

51693205 

0952125 

0 

0 


A(1)=0, 

A(2)=0, 

A(3)=-, 

A(U)=1 

A(5)=0, 

A(6)=0, 

X1  =  0 

X2=FTX 

CALL  NWRPH(X,A,X1,X2) 

Y=XM*X+XB 

RETURN 

END 
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SUBROUTINE  KURVU   TYPE  I 


*ONE  WORD  INTEGERS 
*EXTENDED  PRECISION 
C      R.  A.  BARFORD 

SUBROUTINE  KURVU(X,Y) 
C      UPPER  CURVE  FOR  ACETONE, H20, HCC L3 

Y=((-l.i»99i»i*92*X+0.  91691*23)  *X-  1.1  61*  6217)  *X+0 

RETURN 

END 


99UUI+885 


SUBROUTINE  KURVL   TYPE  I 


*ONE  WORD  INTEGERS 
♦EXTENDED  PRECISION 
C      R.  A.  BARFORD 

SUBROUTINE  KURVL(X,Y) 
C      LOWER  CURVE  FOR  ACETONE, H20, HCCL3 

Y=( ( 1.0952 12 5*X-0. 51693205) *X  +  0. 090264683)  *X  +  0 

RETURN 

END 


005733914 


SUBROUTINE  SBTI E   TYPE  I 


♦ONE  WORD 
♦EXTENDED 


INTEGERS 
PRECISION 


R.  A.  BARFORD 

SUBROUTINE  SBT I E (T I E, NT  I E ) 

DIMENSION  Tl E(20,2) 


TIE(1,1) 
Tl E ( 1, 2  ) 

TIE(2,1) 

TIE(2,2) 

TIE(3,1) 

TIE(3,2) 

TIE(i*,l) 

TIE(U,2) 

TIE(5,1) 

T!E(5,2) 

TIE(6,1) 

TIE(6,2) 

TIE(7,1) 

TIE(7,2) 

Tl £(8,1)  = 

TIE(8,2)  = 

Tl  E(9,l)  =. 

TIE(9,2)  =. 

NTIE  IS  THE 

NTI E  =  9 

RETURN 

END 


2.0 
2.0 

-1.041 
2.1*1*0 
-0.  711*8 
2.306 
-0.5583 
2.608 
-0.1*773 
2.835 
-0.1*586 
3.508 
-0.3858 
1*.  792 
-0.1*350 
II*. 8330 
006 
992 
NUMBER 


OF  TIE  LINES  INCLUDING  THE  BINARY  AXIS  POINTS 
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THE  NEXT  FIVE  SUBROUTINES  DESCRIBE  THE  TYPE  II  SYSTEM  IN  FIG. (2) 


SUBROUTINE  CURVU   TYPE  I  I 


*ONE  WORD 
♦EXTENDED 
C 


A(l) 
A(2) 
A(3) 
A(4) 
A(5) 
A(6) 
XI  = 
X2  = 
CALL 
Y  = 


INTEGERS 
PRECISION 
V.G.  MARTIN 
SUBROUTINE  CURVU ( XM, XB, FTX, X, Y ) 
DIMENSION  A(6) 

=  .95797547  -  XB 

=  -1.2440514  -  XM 

=  1.1377825 

=  -5.0599806 

=  9.3194747 

=-7.3607867 

FTX 

.669 

NWRPH    (X/A/X1/X2) 
XM*X+XB 


RETURN 
END 


SUBROUTINE  CURVL   TYPE  I  I 

♦EXTENDED  PRECISION 
*ONE  WORD  INTEGERS 
C      V.G.  MARTIN 

SUBROUTINE  CURVL ( XM, XB, FTX, X, Y ) 

DIMENSION  A(6) 

A(l)  =  .14126553  -  XB 

A(2)  =  -.55808087  -  XM 

A(3)  =  -1.2103756 

A(4)  =  0.0 

A(5)  =  0.0 

A(6)  =  0.0 

XI  =  0.0 

X2  =  FTX 

CALL  NWRPH(X,A,X1,X2) 

Y  =  XM  *  X  +  XB 

RETURN 

END 
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SUBROUTINE  KURVU   TYPE  I  I 


*ONE  WORD  INTEGERS 
♦EXTENDED  PRECISION 
C      V.G.  MARTIN 

SUBROUTINE  KURVU(X,Y) 

Y  =  ((((-7.3607867  *  X+9 
1-1. 2440514)*X+. 95797547 

RETURN 

END 


3194  74  7)*X-5.0599806)*X+1.1377825)*X 


SUBROUTINE  KURVL   TYPE  I  I 


*ONE  WORD  INTEGERS 
♦EXTENDED  PRECISION 
C      V.G.  MARTIN 

SUBROUTINE  KURVL(X,Y) 


Y  =  (-1 

RETURN 
END 


2103756*X 


55808087)  *X+. 14126553 


SUBROUTINE  SBTIE   TYPE  I  I 


♦ONE  WORD  I NTEGERS 

♦EXTENDED  PRECISION 

C      V.G.  MARTIN 

SUBROUTINE  SBT I E ( T I E, NT  I E ) 
DIMENSION  Tl E(20,2) 


TIE(1,1) 

TIE(1,2) 

TIE(2,1) 

TIE(2,2) 

TIE(3,1) 

TIE(3/2) 

TIE(4,1) 

TIE(4,2) 

TIE(5/1) 

TIE(5/2) 

TIE(6/1) 

TIE(6,2) 

TIE(7,1) 

TIE(7,2) 

TIE(8/1) 

Ti  E ( 8^  2  ) 

NT  I  E  =  8 

RETURN 

END 


.181 

.669 

.0181 

.1104 

.0518 

.4223 

.0594 

1.294 

.0648 

2.675 

.0578 

6.090 

.0678 

38.10 

.14126553 

958 
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MAINLINE  PROGRAM  PUGET 

*ONE  WORD  I NTEGERS 

♦EXTENDED  PRECISION 

*  IOCS (CARD,  TYPEWRITER, KEYBOARD, 11 32  PR  INTER, PAPER  TAPE, DISK) 

**  PUGET  -  VERSION  OF  MAY  20,  1968 

C      V.  G.  MARTIN,  R.  A.  BARFORD 

REAL  LINEA,LINEH 

DIMENSION  WTMEP(5),WTHEX(5),WTACN(5),TI E(20,2) 

DEFINE  FILE  2 ( 2 00, 53, U, NREC ) 

101  FORMAT(IU,2(E13.  5,F7.3),i+E13.5,2F12.3) 

102  FORMATC3F7.U) 

103  FORMATC2X,  'NO.  ',2X,  'WT-CMX  I  N  ' , 2X, ' F -CMX ' , kX, ' WT-CMX  I  N ' , 2X, ' F -CMX 
l',4X,  'WT-CMY  IN1,  4X,'WT-CMY  I  N  ' ,  UX,  '  WT-UPPER  ' ,  5X,  '  WT-  LOWER  ' ,  8X,  '  VO 
2LUME',6X,  'VOLUME'/8X, ' U PPER ' , 5X, ■ U PP ER ' , 5X, 'LOWER', 5 X, ' LOWER',  5X,  ' 
3UPPER",7X,  'LOWER', 3 7X,  ' U PPER ', 7X, ' LOWER  '/  ) 

10U    FORMAT(     'OTOTAL    LOWER    VOLUME    WASHED    OUT    = " ,  F 9  .  k ) 
105    FORMATC     'OTOTAL    CMX    ELUTED    =  '  , F 9 . k ) 

WRITE(3,103) 

READ(2,102)  DENSM,DENSH,DENSA 

TWM  E  P  =  0.0 

TVLOW  =0.0 

NREC  =  1 

CALL  SBTIE(TIE,NTIE) 
1  READ(2 'NREC)WTMEP,WTHEX,WTACN,NDISK 

DO  2  JR  =  1,5 

XX  =  WTMEP(JR)+WTHEX(JR)+WTACN(JR) 

AA  =  WTMEP(JR)/XX  +  .000  5 

BB  =  WTHEX(JR)/XX  +  .0005 

CC  =  WTACN(JR)/XX  +  .  0005 

325  IF  (WTMEP(JR)-(.lE-08)  )  326,  327,  327 

326  WTMEP(JR)  =  0.0 

327  WTTOT  =  WTMEP(JR)  +  WTHEX ( JR ) +WTACN ( JR ) 
FTMEP  =  WTMEP(JR)  /  WTTOT 

FTHEX  =  WTHEX(JR)  /  WTTOT 

FTACN  =  WTACN(JR)  /  WTTOT 

IF(FTMEP)  U0U0,33O,36O 
330  IF(FTHEX-TIE(NTIE,1))520,520,340 
31*0  I  F  CFTHEX-TI  ECNTI  E,2))  350,520,520 
350  FUMEP  =  0.0 

FLMEP  =  0.0 

FUHEX  =  Tl ECNTI E,2) 

FLHEX  =  TIE(NTIE,1) 

GO  TO  500 
370  IF(FTMEP-TIE(1,D)  520,520,375 

360  IF(FTHEX)  ^OUO,  361,  385 

361  IF(TIE(1,1)-1.0)  370,  370,  520 
3  75  IF(FTMEP-TIE(1,2))380,520,520 
380  FUHEX  =  0.0 

FLHEX  =  0.0 
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PUGET    (CONT.) 


PROGRAM 


FUMEP    =    TIE(1,2) 

FLMEP    =    T I E( 1, 1 ) 

GO    TO    500 

CALL    KURVL    (FTMEP,Y) 
C  ANY    POINT    BELOW    Y,    WHICH    IS    ON    THE    LOWER    SOLUBILITY    CURVE,     IS     IN 

C  THE    MISCIBLE    REGION. 

IF(Y-FTHEX)    390,     520,520 

CALL    KURVU(FTMEP,Y) 

IF(FTHEX-Y)    395,520,520 

J    =    1 

J=J  +  1 

TRYY    =    TIE(J,1)    +    Tl E(J,2)*FTMEP 

TEMP    =    FTHEX    -    TRYY 

IF(ABS(TEMP)-.00D    460  ,  40  5,  40  5 

IF  (FTHEX  -  TRYY)  4  10,1160,4  70 

IF  (J-2)  4040,  415,  430 

IF  (Tl E(l,l)-1.0)  420,  4040,  4040 

X  =  -Tl E(J,1)/TI E(J,2) 

GO  TO  440 

L  =  J-l 

TEMP=TIE(L,2)-TI E(J,2) 

X  -   (TIE(J,1)  -  TIE(L,1))/TEMP 

Y  =  Tl E(J,1)  +  TIE(J,2)*X 

TEMP=FTMEP-X 

Tl EM=(FTHEX-Y)/TEMP 

TEMP=TI EM*FTMEP 

Tl EB=FTHEX-TEMP 

CALL    CURVUUI  EM,TIEB,FTMEP,FUMEP,FUHEX) 

CALL    CURVL(TI EM,TI EB, FTMEP, F LMEP, F LHEX ) 

GO    TO    500 

TIEM    =    TIE(J,2) 

Tl EB    =    Tl E(J,1) 

GO    TO    4  50 

I FC J-NT I E+l )    400,480,4040 

X    =    0.0 

GO    TO    440 

FUACN=1.-FUMEP-FUHEX 

FLACN=1.-FLMEP-FLHEX 

TEMP=FUMEP/DENSM+FUHEX/DENSH+FUACN/DENSA 

RHOUP=l./TEMP 

TEMP=FLMEP/DENSM+FLHEX/DENSH+FLACN/DENSA 

RHOLO=l./TEMP 

IF(RHOLO    -    RHOUP)    510,520,530 
510    TEMP    =    FLMEP 

FLMEP    =    FUMEP 

FUMEP    =    TEMP 

TEMP    =    FLHEX 

FLHEX    =    FUHEX 


385 


390 

395 

400 


405 
410 
415 
420 

430 


440 


450 


460 


470 
480 

500 
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PROGRAM    I  I 

PUGET    (CONT.) 

FUHEX    =    TEMP 

TEMP    =    FLACN 

FLACN    =    FUACN 

FUACN    =    TEMP 

TEMP    =    RHOUP 

RHOUP    =    RHOLO 

RHOLO    =    TEMP 

GO    TO    530 
520    FLMEP    =    FTMEP 

FUMEP    =    FLMEP 

FLHEX    =    FTHEX 

FUHEX    =    FLHEX 

FLACN    =    FTACN 

FUACN    =    FLACN 

TEMP=FLMEP/DENSM+FLHEX/DENSH+FLACN/DENSA 

RH0L0=1./TEMP 

WTLOW    =0.0 

FLMEP=0.0 

GO    TO    600 
530    AC    =    FTMEP    -    FLMEP 

BC    =    FTHEX    -    FLHEX 

LINEA    =    SQRT(AC*AC    +    BC*BC) 

AC    =    FUMEP    -    FTMEP 

BC    =    FUHEX    -    FTHEX 

LINEH    =    SQRT(AC*AC    +    BC*BC) 

TEMP    =    LI NEA    +    LI  NEH 

LINEH    =    LINEH    /    TEMP 

WTLOW    =    LI NEH    *    WTTOT 
600    WTUPR    =    WTTOT    -    WTLOW 

WUMEP=WTUPR*FUMEP 

WLMEP=WTLOW*FLMEP 

WUHEX=WTUPR*FUHEX 

WLHEX=WTLOW*FLHEX 

VLOW    =    WTLOW/ RHOLO 

VUP       =    WTUPR/RHOUP 

NFRCT=NDISK-1 

WR I TE ( 3, 10 1 ) NFRCT, WUMEP, FUMEP, WLMEP, F LMEP, WUHEX, WLHEX, WTUPR, WTLOW, 
1 VUP, VLOW 

TVLOW    =    TVLOW    +    VLOW 

TWMEP    =    WTMEP(JR)    +    TWMEP 

2  NDISK  =  NDISK+1 
CALL  DATSW  (3,JSWCH) 
GO  TO  (3,1),JSWCH 

3  WRITE(3,10U)  TVLOW 
WRITE (3, 10 5)  TWMEP 
CALL  EXIT 

1*040  STOP  4040 
END 
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